From c17c8eabc19f942609980b023725baebe42cf259 Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Mon, 3 Nov 2025 09:06:41 +0000 Subject: [PATCH 1/4] bug: split voltage traces before apd90 calculation --- .../AbstractActionPotentialMethod.cpp | 54 +++++++++++-------- 1 file changed, 32 insertions(+), 22 deletions(-) diff --git a/src/single_cell/AbstractActionPotentialMethod.cpp b/src/single_cell/AbstractActionPotentialMethod.cpp index 8b36869..3434f92 100644 --- a/src/single_cell/AbstractActionPotentialMethod.cpp +++ b/src/single_cell/AbstractActionPotentialMethod.cpp @@ -358,15 +358,35 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( // Get voltage properties const unsigned voltage_index = pModel->GetSystemInformation()->GetStateVariableIndex("membrane_voltage"); std::vector voltages = solution.GetVariableAtIndex(voltage_index); - CellProperties voltage_properties(voltages, solution.rGetTimes(), - mActionPotentialThreshold); + std::vector times = solution.rGetTimes(); + std::vector apd90s; std::vector peak_voltages; // See if we can get back some action potential duration(s). try { - apd90s = voltage_properties.GetAllActionPotentialDurations(90); + + // split into num_paces_to_analyze paces for analysis + unsigned current_index = 0; + std::vector voltage_properties; + for (unsigned pace = 0; pace < num_paces_to_analyze; pace++) + { + std::vector pace_voltages; + std::vector pace_times; + const double pace_end_time = (pace + 1) * s1_period; + while (current_index < times.size() && times[current_index] <= pace_end_time) + { + pace_voltages.push_back(voltages[current_index]); + pace_times.push_back(times[current_index]); + current_index++; + } + CellProperties voltage_properties_for_pace(pace_voltages, pace_times, mActionPotentialThreshold); + auto apd90s_for_pace = voltage_properties_for_pace.GetAllActionPotentialDurations(90); + voltage_properties.push_back(voltage_properties_for_pace); + apd90s.insert(std::end(apd90s), std::begin(apd90s_for_pace), std::end(apd90s_for_pace)); + } + if (!mSuppressOutput) { std::cout << "Last " << apd90s.size() @@ -378,25 +398,15 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( std::cout << std::endl; //<< std::flush; } - if (apd90s.size() >= 2u && fabs(apd90s[0] - apd90s[1]) > 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]; - peak_voltages = voltage_properties.GetPeakPotentials(); - rPeak = peak_voltages[0]; - rPeakTime = voltage_properties.GetTimesAtPeakPotentials()[0]; - } - else - { - // Return the last as it is more likely to be the steady state one. - rApd90 = voltage_properties.GetLastActionPotentialDuration(90); - rApd50 = voltage_properties.GetLastActionPotentialDuration(50); - rUpstroke = voltage_properties.GetLastCompleteMaxUpstrokeVelocity(); - rPeak = voltage_properties.GetLastCompletePeakPotential(); - rPeakTime = voltage_properties.GetTimeAtLastCompletePeakPotential(); - } + // if we suspect alternans, analyse the first of the two APs, otherwise the + // second + unsigned analysis_pace_index = (apd90s.size() >= 2u && fabs(apd90s[0] - apd90s[1]) > alternans_threshold) ? 0u : 1u; + rApd90 = voltage_properties[analysis_pace_index].GetLastActionPotentialDuration(90); + rApd50 = voltage_properties[analysis_pace_index].GetLastActionPotentialDuration(50); + rUpstroke = voltage_properties[analysis_pace_index].GetLastCompleteMaxUpstrokeVelocity(); + rPeak = voltage_properties[analysis_pace_index].GetLastCompletePeakPotential(); + rPeakTime = voltage_properties[analysis_pace_index].GetTimeAtLastCompletePeakPotential(); + // It makes sense to return the peak voltage time relative to start of // stimulus application. boost::shared_ptr p_reg_stim = boost::static_pointer_cast( From ad5a334fc32d156f3dd4b1765f3a91e9f2d794bf Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Mon, 3 Nov 2025 09:10:33 +0000 Subject: [PATCH 2/4] tidy up --- src/single_cell/AbstractActionPotentialMethod.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/single_cell/AbstractActionPotentialMethod.cpp b/src/single_cell/AbstractActionPotentialMethod.cpp index 3434f92..dee4a76 100644 --- a/src/single_cell/AbstractActionPotentialMethod.cpp +++ b/src/single_cell/AbstractActionPotentialMethod.cpp @@ -382,8 +382,8 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( current_index++; } CellProperties voltage_properties_for_pace(pace_voltages, pace_times, mActionPotentialThreshold); - auto apd90s_for_pace = voltage_properties_for_pace.GetAllActionPotentialDurations(90); voltage_properties.push_back(voltage_properties_for_pace); + auto apd90s_for_pace = voltage_properties_for_pace.GetAllActionPotentialDurations(90); apd90s.insert(std::end(apd90s), std::begin(apd90s_for_pace), std::end(apd90s_for_pace)); } @@ -400,7 +400,7 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( // if we suspect alternans, analyse the first of the two APs, otherwise the // second - unsigned analysis_pace_index = (apd90s.size() >= 2u && fabs(apd90s[0] - apd90s[1]) > alternans_threshold) ? 0u : 1u; + const unsigned analysis_pace_index = (apd90s.size() >= 2u && fabs(apd90s[0] - apd90s[1]) > alternans_threshold) ? 0u : 1u; rApd90 = voltage_properties[analysis_pace_index].GetLastActionPotentialDuration(90); rApd50 = voltage_properties[analysis_pace_index].GetLastActionPotentialDuration(50); rUpstroke = voltage_properties[analysis_pace_index].GetLastCompleteMaxUpstrokeVelocity(); From 995f22a5bc8eb1941e2cf0b8ca4b87906a37f4ae Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Wed, 5 Nov 2025 16:34:09 +0000 Subject: [PATCH 3/4] use std::find_if and vector contructors --- .../AbstractActionPotentialMethod.cpp | 23 +++++++++++-------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/single_cell/AbstractActionPotentialMethod.cpp b/src/single_cell/AbstractActionPotentialMethod.cpp index dee4a76..aa8257b 100644 --- a/src/single_cell/AbstractActionPotentialMethod.cpp +++ b/src/single_cell/AbstractActionPotentialMethod.cpp @@ -368,23 +368,28 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( { // split into num_paces_to_analyze paces for analysis - unsigned current_index = 0; std::vector voltage_properties; + auto pace_start = std::begin(times); for (unsigned pace = 0; pace < num_paces_to_analyze; pace++) { - std::vector pace_voltages; - std::vector pace_times; + // find the start and end iterators and indices for this pace const double pace_end_time = (pace + 1) * s1_period; - while (current_index < times.size() && times[current_index] <= pace_end_time) - { - pace_voltages.push_back(voltages[current_index]); - pace_times.push_back(times[current_index]); - current_index++; - } + auto pace_end = std::find_if(pace_start, std::end(times), [pace_end_time](double t) { return t > pace_end_time; }); + const size_t pace_start_index = std::distance(std::begin(times), pace_start); + const size_t pace_end_index = std::distance(std::begin(times), pace_end); + + // extract times and voltages + std::vector pace_times(pace_start, pace_end); + std::vector pace_voltages(std::begin(voltages) + pace_start_index, std::begin(voltages) + pace_end_index); + + // contruct CellProperties and save for later CellProperties voltage_properties_for_pace(pace_voltages, pace_times, mActionPotentialThreshold); voltage_properties.push_back(voltage_properties_for_pace); auto apd90s_for_pace = voltage_properties_for_pace.GetAllActionPotentialDurations(90); apd90s.insert(std::end(apd90s), std::begin(apd90s_for_pace), std::end(apd90s_for_pace)); + + // update pace_start for next iteration + pace_start = pace_end; } if (!mSuppressOutput) From f3ab0c1726f52756439cbf87101f8ed2565e4cf2 Mon Sep 17 00:00:00 2001 From: martinjrobins Date: Wed, 5 Nov 2025 16:43:02 +0000 Subject: [PATCH 4/4] we know the pace start and end indices cause its regular --- src/single_cell/AbstractActionPotentialMethod.cpp | 14 ++++---------- 1 file changed, 4 insertions(+), 10 deletions(-) diff --git a/src/single_cell/AbstractActionPotentialMethod.cpp b/src/single_cell/AbstractActionPotentialMethod.cpp index aa8257b..72f417e 100644 --- a/src/single_cell/AbstractActionPotentialMethod.cpp +++ b/src/single_cell/AbstractActionPotentialMethod.cpp @@ -369,17 +369,14 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( // split into num_paces_to_analyze paces for analysis std::vector voltage_properties; - auto pace_start = std::begin(times); for (unsigned pace = 0; pace < num_paces_to_analyze; pace++) { - // find the start and end iterators and indices for this pace - const double pace_end_time = (pace + 1) * s1_period; - auto pace_end = std::find_if(pace_start, std::end(times), [pace_end_time](double t) { return t > pace_end_time; }); - const size_t pace_start_index = std::distance(std::begin(times), pace_start); - const size_t pace_end_index = std::distance(std::begin(times), pace_end); + // find the start and end of the pace + const size_t pace_start_index = static_cast(std::floor(pace * s1_period / printingTimeStep)); + const size_t pace_end_index = static_cast(std::floor((pace + 1) * s1_period / printingTimeStep)); // extract times and voltages - std::vector pace_times(pace_start, pace_end); + std::vector pace_times(std::begin(times) + pace_start_index, std::begin(times) + pace_end_index); std::vector pace_voltages(std::begin(voltages) + pace_start_index, std::begin(voltages) + pace_end_index); // contruct CellProperties and save for later @@ -387,9 +384,6 @@ OdeSolution AbstractActionPotentialMethod::PerformAnalysisOfTwoPaces( voltage_properties.push_back(voltage_properties_for_pace); auto apd90s_for_pace = voltage_properties_for_pace.GetAllActionPotentialDurations(90); apd90s.insert(std::end(apd90s), std::begin(apd90s_for_pace), std::end(apd90s_for_pace)); - - // update pace_start for next iteration - pace_start = pace_end; } if (!mSuppressOutput)