Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
CellProperties.cpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#include <algorithm>
37#include <cassert>
38#include <cmath>
39#include <sstream>
40
42//#include <iostream>
43//#include <iomanip>
44
45#include "CellProperties.hpp"
46#include "Exception.hpp"
47#include "Warnings.hpp"
48
49enum APPhases
50{
51 BELOWTHRESHOLD,
52 ABOVETHRESHOLD
53};
54
56{
57 // Check we have some suitable data to process
58 if (mrTime.size() < 1)
59 {
60 EXCEPTION("Insufficient time steps to calculate physiological properties.");
61 }
62
63 if (mrTime.size() != mrVoltage.size())
64 {
65 EXCEPTION("Time and Voltage series should be the same length. Time.size() = " << mrTime.size() << ", Voltage.size() = " << mrVoltage.size());
66 }
67
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;
77 //boolean to keep track whether we are switching phase from BELOWTHRESHOLD to ABOVETHRESHOLD
78 bool switching_phase = false;
79 bool found_a_flat_bit = false;
80 APPhases ap_phase = BELOWTHRESHOLD;
81
82 unsigned time_steps = mrTime.size() - 1; //The number of time steps is the number of intervals
83
84 double v = mrVoltage[0];
85 double t = mrTime[0];
86 double prev_v = v;
87 double prev_t = t;
88 double voltage_derivative;
89 const double resting_potential_gradient_threshold = 1e-2;
90
91 for (unsigned i = 1; i <= time_steps; i++)
92 {
93 v = mrVoltage[i];
94 t = mrTime[i];
95 voltage_derivative = (t == prev_t) ? 0.0 : (v - prev_v) / (t - prev_t);
96
97 // Look for the max upstroke velocity and when it happens (could be below or above threshold).
98 if (voltage_derivative >= max_upstroke_velocity)
99 {
100 max_upstroke_velocity = voltage_derivative;
101 current_time_of_upstroke_velocity = t;
102 }
103
104 switch (ap_phase)
105 {
106 case BELOWTHRESHOLD:
107 //while below threshold, find the resting value by checking where the velocity is minimal
108 //i.e. when it is flattest. If we can't find a flat bit, instead go for the minimum voltage
109 //seen before the threshold.
110 if (fabs(voltage_derivative) <= current_minimum_velocity && fabs(voltage_derivative) <= resting_potential_gradient_threshold)
111 {
112 current_minimum_velocity = fabs(voltage_derivative);
113 current_resting_value = prev_v;
114 found_a_flat_bit = true;
115 }
116 else if (prev_v < current_resting_value && !found_a_flat_bit)
117 {
118 current_resting_value = prev_v;
119 }
120
121 // If we cross the threshold, this counts as an AP
122 if (v > mThreshold && prev_v <= mThreshold)
123 {
124 //register the resting value and re-initialise the minimum velocity
125 mRestingValues.push_back(current_resting_value);
126 current_minimum_velocity = DBL_MAX;
127 current_resting_value = DBL_MAX;
128
129 //Register the onset time. Linear interpolation.
130 mOnsets.push_back(prev_t + (t - prev_t) / (v - prev_v) * (mThreshold - prev_v));
131
132 //If it is not the first AP, calculate cycle length for the last two APs
133 if (ap_counter > 0)
134 {
135 mCycleLengths.push_back(mOnsets[ap_counter] - mOnsets[ap_counter - 1]);
136 }
137
138 switching_phase = true;
139 found_a_flat_bit = false;
140 ap_phase = ABOVETHRESHOLD;
141 }
142 else
143 {
144 break;
145 }
146 // no break here - deliberate fall through to next case
147
148 case ABOVETHRESHOLD:
149 //While above threshold, look for the peak potential for the current AP
150 if (v > current_peak)
151 {
152 current_peak = v;
153 current_peak_time = t;
154 }
155
156 // we check whether we have above threshold depolarisations
157 // and only if if we haven't just switched from below threshold at this time step.
158 // The latter is to avoid recording things depending on resting behaviour (in case of sudden upstroke from rest)
159 if (prev_voltage_derivative <= 0 && voltage_derivative > 0 && !switching_phase)
160 {
161 counter_of_plateau_depolarisations++;
162 }
163
164 // From the next time step, we are not "switching phase" any longer
165 // (we want to check for above threshold deolarisations)
166 switching_phase = false;
167
168 // If we cross the threshold again, the AP is over
169 // and we register all the parameters.
170 if (v < mThreshold && prev_v >= mThreshold)
171 {
172 // Register peak value for this AP
173 mPeakValues.push_back(current_peak);
174 mTimesAtPeakValues.push_back(current_peak_time);
175 // Re-initialise the current_peak.
176 current_peak = mThreshold;
177 current_peak_time = -DBL_MAX;
178
179 // Register maximum upstroke velocity for this AP
180 mMaxUpstrokeVelocities.push_back(max_upstroke_velocity);
181 // Re-initialise max_upstroke_velocity
182 max_upstroke_velocity = -DBL_MAX;
183
184 // Register time when maximum upstroke velocity occurred for this AP
185 mTimesAtMaxUpstrokeVelocity.push_back(current_time_of_upstroke_velocity);
186 // Re-initialise current_time_of_upstroke_velocity=t;
187 current_time_of_upstroke_velocity = 0.0;
188
189 mCounterOfPlateauDepolarisations.push_back(counter_of_plateau_depolarisations);
190
191 //update the counters.
192 ap_counter++;
193 ap_phase = BELOWTHRESHOLD;
194
195 //reinitialise counter of plateau depolarisations
196 counter_of_plateau_depolarisations = 0;
197 }
198 break;
199 }
200 prev_v = v;
201 prev_t = t;
202 prev_voltage_derivative = voltage_derivative;
203 }
204
205 // One last check. If the simulation ends halfway through an AP
206 // i.e. if the vectors of onsets has more elements than the vectors
207 // of peak and upstroke properties (that are updated at the end of the AP),
208 // then we register the peak and upstroke values so far
209 // for the last incomplete AP.
210 if (mOnsets.size() > mMaxUpstrokeVelocities.size())
211 {
212 mMaxUpstrokeVelocities.push_back(max_upstroke_velocity);
213 mPeakValues.push_back(current_peak);
214 mTimesAtPeakValues.push_back(current_peak_time);
215 mTimesAtMaxUpstrokeVelocity.push_back(current_time_of_upstroke_velocity);
217 }
218}
219
220std::vector<double> CellProperties::CalculateActionPotentialDurations(const double percentage)
221{
223
224 double prev_v = mrVoltage[0];
225 double prev_t = mrTime[0];
226 std::vector<double> apds;
227 std::vector<double> targets;
228 double apd_start_time = DOUBLE_UNSET;
229 double apd_end_time = DOUBLE_UNSET;
230
231 unsigned apd_starting_index = UNSIGNED_UNSET;
232
233 // New algorithm is to loop over APs instead of time.
234 for (unsigned ap_index = 0; ap_index < mPeakValues.size(); ap_index++)
235 {
236 targets.push_back(mRestingValues[ap_index] + 0.01 * (100 - percentage) * (mPeakValues[ap_index] - mRestingValues[ap_index]));
237
238 // We need to look from starting_time_index (just before threshold is reached)
239 unsigned starting_time_index = std::lower_bound(mrTime.begin(), mrTime.end(), mOnsets[ap_index]) - mrTime.begin() - 1u;
240
241 // Now if the target voltage for depolarisation is above the threshold,
242 // we'll need to look forwards in time for the crossing point.
243 if (targets[ap_index] >= mThreshold)
244 {
245 //std::cout << "Looking forwards\n";
246 prev_v = mrVoltage[starting_time_index];
247 prev_t = mrTime[starting_time_index];
248 for (unsigned t = starting_time_index + 1u; t < mrTime.size(); t++)
249 {
250 if (mrVoltage[t] > targets[ap_index])
251 {
252 apd_start_time = prev_t + ((targets[ap_index] - prev_v) / (mrVoltage[t] - prev_v)) * (mrTime[t] - prev_t);
253 apd_starting_index = t;
254 break;
255 }
256 prev_t = mrTime[t];
257 prev_v = mrVoltage[t];
258 }
259 }
260 else // otherwise, we'll need to look backwards.
261 {
262 //std::cout << "Looking backwards\n";
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--)
266 {
267 if (mrVoltage[t] < targets[ap_index])
268 {
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); // Should be a safe conversion since t not allowed to go negative in this loop.
271 break;
272 }
273 prev_t = mrTime[t];
274 prev_v = mrVoltage[t];
275 }
276 }
277
278 // If the start of this AP crossing threshold is before the last one finished,
279 // we are in a wacky regime (see TestCellProperties::TestVeryLongApDetection() for examples of this)
280 // and we need to skip the second one's evaluation.
281 if (apd_end_time != DOUBLE_UNSET && apd_start_time != DOUBLE_UNSET && apd_start_time < apd_end_time)
282 {
283 continue; // Skip to next AP
284 }
285
286 // If there was a previous AP just check for common sense things (to help alert people to above situations)
287 if (ap_index > 0u)
288 {
289 if (fabs(targets[ap_index] - targets[ap_index - 1u]) > 2) // If we see more than a 2mV shift in AP threshold
290 {
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!");
293 }
294 }
295
296 // Now just look forwards for the repolarisation time.
297 if (apd_starting_index != UNSIGNED_UNSET)
298 {
299 prev_t = mrTime[apd_starting_index - 1u];
300 prev_v = mrVoltage[apd_starting_index - 1u];
301 // Now look for a fall below threshold after 1 past the Onset index
302 for (unsigned t = apd_starting_index; t < mrTime.size(); t++)
303 {
304 if (mrVoltage[t] < targets[ap_index])
305 {
306 apd_end_time = mrTime[t - 1] + ((targets[ap_index] - prev_v) / (mrVoltage[t] - prev_v)) * (mrTime[t] - mrTime[t - 1]);
307 apds.push_back(apd_end_time - apd_start_time);
308 break;
309 }
310 prev_t = mrTime[t];
311 prev_v = mrVoltage[t];
312 }
313 }
314 }
315
316 if (apds.size() == 0)
317 {
318 EXCEPTION("No full action potential was recorded");
319 }
320 return apds;
321}
322
324{
326 unsigned size = mPeakValues.size();
327 std::vector<double> amplitudes(size);
328 for (unsigned i = 0; i < size; i++)
329 {
330 amplitudes[i] = (mPeakValues[i] - mRestingValues[i]);
331 }
332 return amplitudes;
333}
334
339
341{
342 // mOnsets and mRestingValues are all
343 // set at the same time, so checking one should suffice.
344 if (mOnsets.empty())
345 {
346 // possible false error here if the simulation started at time < 0
347 EXCEPTION("AP did not occur, never exceeded threshold voltage.");
348 }
349}
350
352{
353 // mPeakValues, mMaxUpstrokeVelocities and mTimesAtMaxUpstrokeVelocity are all
354 // set at the same time so checking one should suffice.
355 if (mMaxUpstrokeVelocities.empty())
356 {
357 EXCEPTION("AP did not occur, never descended past threshold voltage.");
358 }
359}
360
361//
362// The Get std::vector<double> methods
363//
364
366{
367 return mCycleLengths;
368}
369
371{
373 return mRestingValues;
374}
375
377{
379 return mPeakValues;
380}
381
387
393
399
400std::vector<double> CellProperties::GetAllActionPotentialDurations(const double percentage)
401{
402 return CalculateActionPotentialDurations(percentage);
403}
404
409//
410// The Get <double> methods
411//
412
418
420{
422 double peak_value;
423 if (mUnfinishedActionPotentials && mPeakValues.size() > 1u)
424 {
425 peak_value = mPeakValues[mPeakValues.size() - 2u];
426 }
427 else if (mUnfinishedActionPotentials && mPeakValues.size() == 1u)
428 {
429 EXCEPTION("No peak potential matching a full action potential was recorded.");
430 }
431 else
432 {
433 peak_value = mPeakValues.back();
434 }
435 return peak_value;
436}
437
439{
441 double peak_value_time;
443 {
444 peak_value_time = mTimesAtPeakValues[mTimesAtPeakValues.size() - 2u];
445 }
446 else if (mUnfinishedActionPotentials && mTimesAtPeakValues.size() == 1u)
447 {
448 EXCEPTION("No peak potential matching a full action potential was recorded.");
449 }
450 else
451 {
452 peak_value_time = mTimesAtPeakValues.back();
453 }
454 return peak_value_time;
455}
456
462
464{
466 double max_upstroke;
468 {
469 max_upstroke = mMaxUpstrokeVelocities[mMaxUpstrokeVelocities.size() - 2u];
470 }
472 {
473 EXCEPTION("No MaxUpstrokeVelocity matching a full action potential was recorded.");
474 }
475 else
476 {
477 max_upstroke = mMaxUpstrokeVelocities.back();
478 }
479 return max_upstroke;
480}
481
487
493
495{
497 double max_upstroke_time;
499 {
500 max_upstroke_time = mTimesAtMaxUpstrokeVelocity[mTimesAtMaxUpstrokeVelocity.size() - 2u];
501 }
503 {
504 EXCEPTION("No TimeAtMaxUpstrokeVelocity matching a full action potential was recorded.");
505 }
506 else
507 {
508 max_upstroke_time = mTimesAtMaxUpstrokeVelocity.back();
509 }
510 return max_upstroke_time;
511}
512
514{
515 // We tested this returns a non-empty vector in the method.
516 return CalculateActionPotentialDurations(percentage).back();
517}
518
523
525{
526 // We tested something sensible has happened in the vector returning method.
527 return GetRestingPotentials().back();
528}
const double DOUBLE_UNSET
Definition Exception.hpp:57
#define EXCEPTION(message)
const unsigned UNSIGNED_UNSET
Definition Exception.hpp:53
std::vector< double > mOnsets
std::vector< double > GetTimesAtMaxUpstrokeVelocity()
double GetTimeAtLastPeakPotential()
unsigned GetNumberOfAboveThresholdDepolarisationsForLastAp()
std::vector< double > GetPeakPotentials()
double GetLastRestingPotential()
const std::vector< double > & mrTime
std::vector< double > GetCycleLengths()
std::vector< double > mPeakValues
void CheckExceededThreshold(void)
double GetTimeAtLastCompleteMaxUpstrokeVelocity()
std::vector< unsigned > GetNumberOfAboveThresholdDepolarisationsForAllAps()
std::vector< double > mTimesAtPeakValues
std::vector< double > mMaxUpstrokeVelocities
double GetLastActionPotentialAmplitude()
std::vector< unsigned > mCounterOfPlateauDepolarisations
std::vector< double > GetAllActionPotentialDurations(const double percentage)
std::vector< double > mRestingValues
double GetTimeAtLastCompletePeakPotential()
std::vector< double > CalculateActionPotentialDurations(const double percentage)
void CheckReturnedToThreshold(void)
std::vector< double > mTimesAtMaxUpstrokeVelocity
double GetLastPeakPotential()
std::vector< double > GetTimesAtPeakPotentials()
double GetLastCompleteMaxUpstrokeVelocity()
std::vector< double > GetMaxUpstrokeVelocities()
bool mUnfinishedActionPotentials
std::vector< double > GetActionPotentialAmplitudes()
const std::vector< double > & mrVoltage
std::vector< double > mCycleLengths
double GetTimeAtLastMaxUpstrokeVelocity()
double GetLastMaxUpstrokeVelocity()
double GetLastActionPotentialDuration(const double percentage)
double GetLastCompletePeakPotential()
std::vector< double > GetRestingPotentials()