Chaste  Release::2018.1
CellProperties.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF 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 
49 enum 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 
220 std::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 
336 {
337  return GetActionPotentialAmplitudes().back();
338 }
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 
365 std::vector<double> CellProperties::GetCycleLengths()
366 {
367  return mCycleLengths;
368 }
369 
371 {
373  return mRestingValues;
374 }
375 
377 {
379  return mPeakValues;
380 }
381 
383 {
385  return mTimesAtPeakValues;
386 }
387 
389 {
391  return mMaxUpstrokeVelocities;
392 }
393 
395 {
398 }
399 
400 std::vector<double> CellProperties::GetAllActionPotentialDurations(const double percentage)
401 {
402  return CalculateActionPotentialDurations(percentage);
403 }
404 
406 {
408 }
409 //
410 // The Get <double> methods
411 //
412 
414 {
416  return mPeakValues.back();
417 }
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 
458 {
460  return mMaxUpstrokeVelocities.back();
461 }
462 
464 {
466  double max_upstroke;
468  {
469  max_upstroke = mMaxUpstrokeVelocities[mMaxUpstrokeVelocities.size() - 2u];
470  }
471  else if (mUnfinishedActionPotentials && mMaxUpstrokeVelocities.size() == 1u)
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 
483 {
485  return mTimesAtMaxUpstrokeVelocity.back();
486 }
487 
489 {
491  return mTimesAtPeakValues.back();
492 }
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 
513 double CellProperties::GetLastActionPotentialDuration(const double percentage)
514 {
515  // We tested this returns a non-empty vector in the method.
516  return CalculateActionPotentialDurations(percentage).back();
517 }
518 
520 {
521  return mCounterOfPlateauDepolarisations.back();
522 }
523 
525 {
526  // We tested something sensible has happened in the vector returning method.
527  return GetRestingPotentials().back();
528 }
std::vector< double > GetAllActionPotentialDurations(const double percentage)
double GetLastMaxUpstrokeVelocity()
std::vector< double > mTimesAtMaxUpstrokeVelocity
double GetTimeAtLastCompletePeakPotential()
double GetTimeAtLastCompleteMaxUpstrokeVelocity()
std::vector< double > GetRestingPotentials()
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::vector< double > GetMaxUpstrokeVelocities()
std::vector< double > GetTimesAtPeakPotentials()
const std::vector< double > & mrTime
double GetLastPeakPotential()
double GetTimeAtLastPeakPotential()
std::vector< double > GetActionPotentialAmplitudes()
std::vector< double > CalculateActionPotentialDurations(const double percentage)
unsigned GetNumberOfAboveThresholdDepolarisationsForLastAp()
const double DOUBLE_UNSET
Definition: Exception.hpp:56
std::vector< double > mCycleLengths
bool mUnfinishedActionPotentials
std::vector< double > GetTimesAtMaxUpstrokeVelocity()
double GetLastActionPotentialDuration(const double percentage)
double GetLastCompleteMaxUpstrokeVelocity()
double GetLastActionPotentialAmplitude()
const unsigned UNSIGNED_UNSET
Definition: Exception.hpp:52
void CheckExceededThreshold(void)
std::vector< double > mPeakValues
std::vector< double > mMaxUpstrokeVelocities
std::vector< unsigned > GetNumberOfAboveThresholdDepolarisationsForAllAps()
std::vector< double > mOnsets
std::vector< unsigned > mCounterOfPlateauDepolarisations
std::vector< double > mTimesAtPeakValues
std::vector< double > GetCycleLengths()
double GetLastCompletePeakPotential()
double GetLastRestingPotential()
void CalculateProperties()
std::vector< double > mRestingValues
std::vector< double > GetPeakPotentials()
void CheckReturnedToThreshold(void)
const std::vector< double > & mrVoltage
double GetTimeAtLastMaxUpstrokeVelocity()