diff --git a/src/single_cell/AbstractActionPotentialMethod.cpp b/src/single_cell/AbstractActionPotentialMethod.cpp index 8b36869..72f417e 100644 --- a/src/single_cell/AbstractActionPotentialMethod.cpp +++ b/src/single_cell/AbstractActionPotentialMethod.cpp @@ -358,15 +358,34 @@ 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 + std::vector voltage_properties; + for (unsigned pace = 0; pace < num_paces_to_analyze; pace++) + { + // 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(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 + 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)); + } + if (!mSuppressOutput) { std::cout << "Last " << apd90s.size() @@ -378,25 +397,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 + 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(); + 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(