60 EXCEPTION(
"Insufficient time steps to calculate physiological properties.");
65 EXCEPTION(
"Time and Voltage series should be the same length. Time.size() = " <<
mrTime.size() <<
", Voltage.size() = " <<
mrVoltage.size());
68 double max_upstroke_velocity = -DBL_MAX;
69 double current_time_of_upstroke_velocity = 0;
70 double current_resting_value = DBL_MAX;
71 double current_peak = -DBL_MAX;
72 double current_peak_time = -DBL_MAX;
73 double current_minimum_velocity = DBL_MAX;
74 double prev_voltage_derivative = 0;
75 unsigned ap_counter = 0;
76 unsigned counter_of_plateau_depolarisations = 0;
78 bool switching_phase =
false;
79 bool found_a_flat_bit =
false;
80 APPhases ap_phase = BELOWTHRESHOLD;
82 unsigned time_steps =
mrTime.size() - 1;
88 double voltage_derivative;
89 const double resting_potential_gradient_threshold = 1e-2;
91 for (
unsigned i = 1; i <= time_steps; i++)
95 voltage_derivative = (t == prev_t) ? 0.0 : (v - prev_v) / (t - prev_t);
98 if (voltage_derivative >= max_upstroke_velocity)
100 max_upstroke_velocity = voltage_derivative;
101 current_time_of_upstroke_velocity = t;
110 if (fabs(voltage_derivative) <= current_minimum_velocity && fabs(voltage_derivative) <= resting_potential_gradient_threshold)
112 current_minimum_velocity = fabs(voltage_derivative);
113 current_resting_value = prev_v;
114 found_a_flat_bit =
true;
116 else if (prev_v < current_resting_value && !found_a_flat_bit)
118 current_resting_value = prev_v;
126 current_minimum_velocity = DBL_MAX;
127 current_resting_value = DBL_MAX;
138 switching_phase =
true;
139 found_a_flat_bit =
false;
140 ap_phase = ABOVETHRESHOLD;
150 if (v > current_peak)
153 current_peak_time = t;
159 if (prev_voltage_derivative <= 0 && voltage_derivative > 0 && !switching_phase)
161 counter_of_plateau_depolarisations++;
166 switching_phase =
false;
177 current_peak_time = -DBL_MAX;
182 max_upstroke_velocity = -DBL_MAX;
187 current_time_of_upstroke_velocity = 0.0;
193 ap_phase = BELOWTHRESHOLD;
196 counter_of_plateau_depolarisations = 0;
202 prev_voltage_derivative = voltage_derivative;
225 double prev_t =
mrTime[0];
226 std::vector<double> apds;
227 std::vector<double> targets;
234 for (
unsigned ap_index = 0; ap_index <
mPeakValues.size(); ap_index++)
247 prev_t =
mrTime[starting_time_index];
248 for (
unsigned t = starting_time_index + 1u; t <
mrTime.size(); t++)
252 apd_start_time = prev_t + ((targets[ap_index] - prev_v) / (
mrVoltage[t] - prev_v)) * (
mrTime[t] - prev_t);
253 apd_starting_index = t;
263 prev_v =
mrVoltage[starting_time_index + 1];
264 prev_t =
mrTime[starting_time_index + 1];
265 for (
int t = starting_time_index; t >= 0; t--)
269 apd_start_time = prev_t + ((targets[ap_index] - prev_v) / (
mrVoltage[t] - prev_v)) * (
mrTime[t] - prev_t);
270 apd_starting_index = (
unsigned)(t + 1);
289 if (fabs(targets[ap_index] - targets[ap_index - 1u]) > 2)
291 WARNING(
"The voltage threshold for measuring APD" << percentage <<
" changed from "
292 << targets[ap_index - 1u] <<
" to " << targets[ap_index] <<
" in this AP trace, which may suggest CellProperties isn't telling you anything very sensible!");
299 prev_t =
mrTime[apd_starting_index - 1u];
300 prev_v =
mrVoltage[apd_starting_index - 1u];
302 for (
unsigned t = apd_starting_index; t <
mrTime.size(); t++)
307 apds.push_back(apd_end_time - apd_start_time);
316 if (apds.size() == 0)
318 EXCEPTION(
"No full action potential was recorded");