feat: run 3 paces to calculate apd90

CMakeLists.txt

@@ -42,7 +42,7 @@
 #
 # Note that the order in which components are specified does not matter.
 
-set(ApPredict_CELLML_TAG e394aba8a7ca50f72ee212e21431045eae8b5ca4)
+set(ApPredict_CELLML_TAG fb76f77ef973ff4e18c2b36b841ff3be5e59b5c8)
 set(ApPredict_CELLML_URL "https://raw.githubusercontent.com/Chaste/cellml/${ApPredict_CELLML_TAG}/cellml")
 set(ApPredict_CELLML_SRC "${CMAKE_CURRENT_SOURCE_DIR}/src/cellml")
 set(ApPredict_CELLML_FILES
@@ -59,6 +59,7 @@ set(ApPredict_CELLML_FILES
     priebe_beuckelmann_1998.cellml
     shannon_wang_puglisi_weber_bers_2004.cellml
     ten_tusscher_model_2006_epi.cellml
+    noble_model_1962.cellml
 )
 
 message (STATUS "Downloading ApPredict cellml files to ${ApPredict_CELLML_SRC}")

src/fortests/ActionPotentialDownsampler.cpp

@@ -46,7 +46,8 @@ ActionPotentialDownsampler::ActionPotentialDownsampler(const std::string& rFolde
                                                        const std::vector<double>& rTimes,
                                                        const std::vector<double>& rVoltages,
                                                        double window,
-                                                       double stimulusStart)
+                                                       double stimulusStart,
+                                                       double windowStart)
 {
     OutputFileHandler handler(rFoldername, false); // don't wipe the folder!
     out_stream output_file = handler.OpenOutputFile(rFilename);
@@ -59,8 +60,10 @@ ActionPotentialDownsampler::ActionPotentialDownsampler(const std::string& rFolde
 
     assert(rTimes.size()>0);
     double start_time_for_this_pace = rTimes[0] + stimulusStart;
+    double dt = rTimes[1] - rTimes[0];
+    unsigned start_index = static_cast<unsigned>(floor(windowStart / dt));
 
-    for (unsigned i=0; i<rVoltages.size(); i++)
+    for (unsigned i=start_index; i<rVoltages.size(); i++)
     {
         if (rTimes[i] - start_time_for_this_pace > window)
         {
@@ -71,7 +74,7 @@ ActionPotentialDownsampler::ActionPotentialDownsampler(const std::string& rFolde
         {
             // A new bit of code to downsample the output so flot doesn't timeout loading large datasets
             // We want to plot the first point, but not necessarily the last if we are repolarized.
-            if (i>0 && (i<rVoltages.size() || rVoltages[i] < -50)) // if we aren't at the beginning or the end of the trace
+            if (i>start_index && (i<rVoltages.size() || rVoltages[i] < -50)) // if we aren't at the beginning or the end of the trace
             {
             	if (fabs(rVoltages[i] - last_voltage_printed) > 1.0 /*mV*/ )
                 {

src/fortests/ActionPotentialDownsampler.hpp

@@ -62,7 +62,8 @@ class ActionPotentialDownsampler
                                const std::vector<double>& rTimes,
                                const std::vector<double>& rVoltages,
                                double window,
-                               double stimulusStart=0);
+                               double stimulusStart=0,
+                               double windowStart=0);
 };
 
 #endif // ACTIONPOTENTIALDOWNSAMPLER_HPP_

src/fortests/SetupModel.cpp

@@ -62,7 +62,8 @@ const std::unordered_set<std::string> SetupModel::forceNumericalJModels = { "hun
 
 SetupModel::SetupModel(const double& rHertz,
                        unsigned modelIndex,
-                       boost::shared_ptr<OutputFileHandler> pHandler)
+                       boost::shared_ptr<OutputFileHandler> pHandler,
+                       bool useSuggestedCycleLength)
         : mpHandler(pHandler)
 {
     // Cvode cells use a CVODE solver regardless of which standard solver is passed in.
@@ -149,7 +150,7 @@ SetupModel::SetupModel(const double& rHertz,
         s_magnitude = p_reg_stim->GetMagnitude();
         s_duration = p_reg_stim->GetDuration();
     }
-    else if (mpModel->HasAttribute("SuggestedCycleLength"))
+    else if (mpModel->HasAttribute("SuggestedCycleLength") && useSuggestedCycleLength)
     {
         // If the model has no stimulus current it can be given this attribute tag
         // so we know roughly what cycle length to use on it...
@@ -195,4 +196,4 @@ SetupModel::SetupModel(const double& rHertz,
 boost::shared_ptr<AbstractCvodeCell> SetupModel::GetModel()
 {
     return mpModel;
-}
+}

src/fortests/SetupModel.hpp

@@ -91,7 +91,8 @@ class SetupModel
      */
     SetupModel(const double& rHertz,
                unsigned model_index = UNSIGNED_UNSET,
-               boost::shared_ptr<OutputFileHandler> pHandler = boost::shared_ptr<OutputFileHandler>());
+               boost::shared_ptr<OutputFileHandler> pHandler = boost::shared_ptr<OutputFileHandler>(),
+               bool useSuggestedCycleLength = true);
 
     /**
      * Print the option and possible arguments that this class takes.

src/single_cell/AbstractActionPotentialMethod.cpp

@@ -306,9 +306,6 @@ OdeSolution AbstractActionPotentialMethod::SteadyStatePacingExperiment(
     if (mRepeat)
     {
         // std::cout << "Repeating simulation to order alternans APs consistently...\n";
-        // If we might benefit from pushing forward one period and re-analysing...
-        PushModelForwardOneS1Interval(pModel, s1_period, maximum_time_step);
-
         solution = PerformAnalysisOfTwoPaces(
             pModel, rApd90, rApd50, rUpstroke, rPeak, rPeakTime, rCaMax, rCaMin,
             s1_period, maximum_time_step, printingTimeStep, conc);
@@ -341,7 +338,7 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
     // We usually analyse two paces to look for alternans.
     // Have observed three period behaviour or more, so this is a hardcoded option
     // for now.
-    const unsigned num_paces_to_analyze = 2u;
+    const unsigned num_paces_to_analyze = 3u;
     mRepeat = false;
     const double alternans_threshold = 1; // ms in APD90 - hardcoded,
     // could make an option in future. But if it is any smaller alternans 'comes
@@ -378,15 +375,15 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
             std::cout << std::endl; //<< std::flush;
         }
 
-        if (apd90s.size() >= 2u && fabs(apd90s[0] - apd90s[1]) > alternans_threshold)
+        if (apd90s.size() >= 3u && fabs(apd90s[1] - apd90s[2]) > alternans_threshold)
         {
             // We suspect alternans, and analyse the first of the two APs
-            rApd90 = voltage_properties.GetAllActionPotentialDurations(90)[0];
-            rApd50 = voltage_properties.GetAllActionPotentialDurations(50)[0];
-            rUpstroke = voltage_properties.GetMaxUpstrokeVelocities()[0];
+            rApd90 = voltage_properties.GetAllActionPotentialDurations(90)[1];
+            rApd50 = voltage_properties.GetAllActionPotentialDurations(50)[1];
+            rUpstroke = voltage_properties.GetMaxUpstrokeVelocities()[1];
             peak_voltages = voltage_properties.GetPeakPotentials();
-            rPeak = peak_voltages[0];
-            rPeakTime = voltage_properties.GetTimesAtPeakPotentials()[0];
+            rPeak = peak_voltages[1];
+            rPeakTime = voltage_properties.GetTimesAtPeakPotentials()[1];
         }
         else
         {
@@ -480,10 +477,10 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
         assert(apd90s.size() > 0u);
 
         // Deal with the case when we are in alternans
-        if (apd90s.size() >= 2u && fabs(apd90s[0] - apd90s[1]) > alternans_threshold)
+        if (apd90s.size() >= 3u && fabs(apd90s[1] - apd90s[2]) > alternans_threshold)
         {
             mPeriodTwoBehaviour = true;
-            if (apd90s[1] > apd90s[0] && mRepeatNumber == 0u)
+            if (apd90s[2] > apd90s[1] && mRepeatNumber == 0u)
             {
                 // Redo so that we always plot the longest AP first.
                 mRepeat = true;
@@ -500,7 +497,7 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
                     // 2. That this isn't caused by a slow depolarisation which tends to
                     // give varied Peak Vm
                     // (see test case 9 in O'Hara model TestTroublesomeApEvaluations)
-                    if (mDefaultParametersApd90 != DOUBLE_UNSET && (apd90s[0] > mDefaultParametersApd90 && apd90s[1] > mDefaultParametersApd90) && (fabs(peak_voltages[0] - peak_voltages[1]) < 10 /*mV*/)) // magic number!
+                    if (mDefaultParametersApd90 != DOUBLE_UNSET && (apd90s[1] > mDefaultParametersApd90 && apd90s[2] > mDefaultParametersApd90) && (fabs(peak_voltages[1] - peak_voltages[2]) < 10 /*mV*/)) // magic number!
                     {
                         // We have alternans tending to long/no repolarisation
                         mErrorCode = 6u;
@@ -519,8 +516,8 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
                 {
                     message << "At a concentration of " << conc << "uM: ";
                 }
-                message << "possible alternans detected, APD90s = " << apd90s[0] << ", "
-                        << apd90s[1] << " ms";
+                message << "possible alternans detected, APD90s = " << apd90s[1] << ", "
+                        << apd90s[2] << " ms";
                 std::string message_string = message.str();
                 WriteMessageToFile(message_string);
             }
@@ -535,16 +532,20 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
             }
         }
 
-        // Deal with the case when there was only one AP but two stimuli
+        // Deal with the case when there was only less than three AP but three stimuli
         if (apd90s.size() < num_paces_to_analyze && mSuccessful)
         {
             std::stringstream message;
             if (conc != DOUBLE_UNSET)
             {
                 message << "At a concentration of " << conc << "uM: ";
             }
-            message << "only one action potential was recorded (" << apd90s[0]
-                    << " ms) for two paces of " << s1_period << "ms.";
+            message << "less than three action potential was recorded (";
+            for (unsigned i = 0; i < apd90s.size(); i++)
+            {
+                message << apd90s[i] << " ms, ";
+            }
+            message << ") for three paces of " << s1_period << "ms.";
             std::string message_string = message.str();
             WriteMessageToFile(message_string);
 
@@ -555,7 +556,15 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces(
                 // absence of this caused a bit of a bug where a period 3
                 // non-repolarising looked like non
                 // depolarising...
-                if (apd90s[0] > s1_period || voltages.back() >= mActionPotentialThreshold)
+                bool an_ap_greater_than_period = false;
+                for (unsigned i = 0; i < apd90s.size(); i++)
+                {
+                    if (apd90s[i] > s1_period)
+                    {
+                        an_ap_greater_than_period = true;
+                    }
+                }
+                if (an_ap_greater_than_period || voltages.back() >= mActionPotentialThreshold)
                 {
                     mErrorCode = 3u;
                     mErrorMessage = "NoActionPotential_3";

src/single_cell/ApPredictMethods.cpp

@@ -1520,15 +1520,16 @@ void ApPredictMethods::CommonRunMethod()
         boost::shared_ptr<RegularStimulus> p_default_stimulus = boost::static_pointer_cast<RegularStimulus>(
             mpModel->GetStimulusFunction());
         double s1_period = p_default_stimulus->GetPeriod();
-        double s_start = p_default_stimulus->GetStartTime();
+        double s_start = p_default_stimulus->GetStartTime() + s1_period;
         std::vector<double> voltages = solution.GetVariableAtIndex(mpModel->GetSystemInformation()->GetStateVariableIndex("membrane_voltage"));
         double window = s1_period;
         if (this->mPeriodTwoBehaviour)
         {
             window *= 2.0;
         }
+        double data_start = s1_period;
         ActionPotentialDownsampler(mOutputFolder, filename.str(),
-                                   solution.rGetTimes(), voltages, window, s_start);
+                                   solution.rGetTimes(), voltages, window, s_start, data_start);
     } // Conc
 
     if (!reliable_credible_intervals)

test/TestApPredict.hpp

@@ -127,8 +127,8 @@ class TestApPredict : public CxxTest::TestSuite
 
     void TestVoltageThresholdDetectionAlgorithm()
     {
-        std::vector<double> thresholds_for_each_model = boost::assign::list_of(-46.7750) /*Shannon etc.*/
-            (-23.0772)(-34.6525)(-35.9230)(-28.4091)(-38.4384)(-39.1051)(-38.636)(-26.9059);
+        std::vector<double> thresholds_for_each_model = boost::assign::list_of(-46.7583) /*Shannon etc.*/
+            (-23.0772)(-34.6375)(-35.9230)(-28.4091)(-38.4384)(-39.1206)(-38.636)(-26.9059);
 
         for (unsigned model_index = 1; model_index < 10u; model_index++)
         {
@@ -206,7 +206,7 @@ class TestApPredict : public CxxTest::TestSuite
             std::vector<double> apd90s = methods.GetApd90s();
             TS_ASSERT_EQUALS(apd90s.size(), 2u);
             TS_ASSERT_DELTA(apd90s[0], 198.381, 2e-2);
-            TS_ASSERT_DELTA(apd90s[1], 198.3093, 2e-2);
+            TS_ASSERT_DELTA(apd90s[1], 198.2879, 2e-2);
         }
     }
 

test/TestTroublesomeApEvaluations.hpp

@@ -39,6 +39,7 @@ OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #include <cxxtest/TestSuite.h>
 #include "SetupModel.hpp"
 #include "SingleActionPotentialPrediction.hpp"
+#include "CommandLineArgumentsMocker.hpp"
 #include "ZeroStimulus.hpp"
 
 class TestTroublesomeApEvaluations : public CxxTest::TestSuite
@@ -94,6 +95,16 @@ class TestTroublesomeApEvaluations : public CxxTest::TestSuite
     }
 
 public:
+    void TestTroublesomeNobleCases()
+    {
+        CommandLineArgumentsMocker wrapper("--model noble_model_1962 --pacing-stim-magnitude -25.5");
+        SetupModel setup(2.0, UNSIGNED_UNSET, boost::shared_ptr<OutputFileHandler>(), false); // Noble 1962, 2.0Hz
+        boost::shared_ptr<AbstractCvodeCell> p_model = setup.GetModel();
+        std::string message = Run(p_model, 100, true, DOUBLE_UNSET, DOUBLE_UNSET, "_noble_normal");
+        // Noble 1962 model has a non-flat resting potential, so we just check that no error occurs in the AP evaluation.
+        TS_ASSERT_EQUALS(message, "No error");
+    }
+
     void TestTroublesomeTenTusscherCases()
     {
         SetupModel setup(1.0, 2u); // TT06, 1.0Hz