CellProperties.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include <cmath>
00030 #include <sstream>
00031 #include <cassert>
00032 
00033 #include "CellProperties.hpp"
00034 #include "Exception.hpp"
00035 
00036 
00037 enum APPhases { BELOWTHRESHOLD , ABOVETHRESHOLD };
00038 
00039 void CellProperties::CalculateProperties()
00040 {
00041     // Check we have some suitable data to process
00042     if (mrTime.size() < 1)
00043     {
00044         EXCEPTION("Insufficient time steps to calculate physiological properties.");
00045     }
00046 
00047     if (mrTime.size() != mrVoltage.size())
00048     {
00049         std::stringstream exception_message;
00050         exception_message << "Time and Voltage series should be the same length. Time.size() = " << mrTime.size() << ", Voltage.size() = " << mrVoltage.size();
00051         EXCEPTION(exception_message.str());
00052     }
00053 
00054 
00055     double max_upstroke_velocity = -DBL_MAX;
00056     double current_time_of_upstroke_velocity = 0;
00057     double current_resting_value=DBL_MAX;
00058     double current_peak=-DBL_MAX;
00059     double current_minimum_velocity=DBL_MAX;
00060     double prev_voltage_derivative=0;
00061     unsigned ap_counter = 0;
00062     unsigned counter_of_plateau_depolarisations = 0;
00063     //boolean to keep track whether we are switching phase from BELOWTHRESHOLD to ABOVETHRESHOLD
00064     bool switching_phase = false;
00065     bool found_a_flat_bit=false;
00066     APPhases ap_phase = BELOWTHRESHOLD;
00067 
00068     unsigned time_steps = mrTime.size()-1; //The number of time steps is the number of intervals
00069 
00070     double v = mrVoltage[0];
00071     double t = mrTime[0];
00072     double prev_v = v;
00073     double prev_t = t;
00074     double voltage_derivative;
00075     const double resting_potential_gradient_threshold = 1e-2; 
00076 
00077     for (unsigned i=1; i<=time_steps; i++)
00078     {
00079         v = mrVoltage[i];
00080         t = mrTime[i];
00081         voltage_derivative = (t == prev_t) ? 0.0 : (v - prev_v) / (t - prev_t);
00082 
00083         // Look for the max upstroke velocity and when it happens (could be below or above threshold).
00084         if (voltage_derivative >= max_upstroke_velocity)
00085         {
00086             max_upstroke_velocity = voltage_derivative;
00087             current_time_of_upstroke_velocity = t;
00088 
00089         }
00090 
00091         switch (ap_phase)
00092         {
00093             case BELOWTHRESHOLD:
00094                 //while below threshold, find the resting value by checking where the velocity is minimal
00095                 //i.e. when it is flattest. If we can't find a flat bit, instead go for the minimum voltage
00096                 //seen before the threshold.
00097                 if (fabs(voltage_derivative)<=current_minimum_velocity && fabs(voltage_derivative)<=resting_potential_gradient_threshold)
00098                 {
00099                     current_minimum_velocity=fabs(voltage_derivative);
00100                     current_resting_value = prev_v;
00101                     found_a_flat_bit=true;
00102                 }
00103                 else if(prev_v < current_resting_value && !found_a_flat_bit)
00104                 {
00105                     current_resting_value = prev_v;
00106                 }
00107 
00108                 // If we cross the threshold, this counts as an AP
00109                 if ( v>mThreshold && prev_v <= mThreshold )
00110                 {
00111                     //register the resting value and re-initialise the minimum velocity
00112                     mRestingValues.push_back(current_resting_value);
00113                     current_minimum_velocity = DBL_MAX;
00114                     current_resting_value = DBL_MAX;
00115 
00116                     //Register the onset time. Linear interpolation.
00117                     mOnsets.push_back(prev_t + (t-prev_t)/(v-prev_v)*(mThreshold-prev_v));
00118 
00119                     //If it is not the first AP, calculate cycle length for the last two APs
00120                     if (ap_counter>0)
00121                     {
00122                         mCycleLengths.push_back( mOnsets[ap_counter]-mOnsets[ap_counter-1] );
00123                     }
00124 
00125                     switching_phase = true;
00126                     found_a_flat_bit = false;
00127                     ap_phase = ABOVETHRESHOLD;
00128                     // no break here - deliberate fall through to next case
00129                 }
00130                 else
00131                 {
00132                     break;
00133                 }
00134 
00135             case ABOVETHRESHOLD:
00136                 //While above threshold, look for the peak potential for the current AP
00137                 if (v>current_peak)
00138                 {
00139                    current_peak = v;
00140                 }
00141 
00142                 // we check whether we have above threshold depolarisations
00143                 // and only if if we haven't just switched from below threshold at this time step.
00144                 // The latter is to avoid recording things depending on resting behaviour (in case of sudden upstroke from rest)
00145                 if (prev_voltage_derivative<=0 && voltage_derivative>0 && !switching_phase)
00146                 {
00147                     counter_of_plateau_depolarisations++;
00148                 }
00149 
00150                 // From the next time step, we are not "switching phase" any longer
00151                 // (we want to check for above threshold deolarisations)
00152                 switching_phase = false;
00153 
00154                 // If we cross the threshold again, the AP is over
00155                 // and we register all the parameters.
00156                 if ( v<mThreshold && prev_v >= mThreshold )
00157                 {
00158                     //register peak value for this AP
00159                     mPeakValues.push_back(current_peak);
00160                     //Re-initialise the current_peak.
00161                     current_peak = mThreshold;
00162 
00163                     //register maximum upstroke velocity for this AP
00164                     mMaxUpstrokeVelocities.push_back(max_upstroke_velocity);
00165                     //re-initialise max_upstroke_velocity
00166                     max_upstroke_velocity = -DBL_MAX;
00167 
00168                     //register time when maximum upstroke velocity occurred for this AP
00169                     mTimesAtMaxUpstrokeVelocity.push_back(current_time_of_upstroke_velocity);
00170                     //re-initialise current_time_of_upstroke_velocity=t;
00171                     current_time_of_upstroke_velocity = 0.0;
00172 
00173                     mCounterOfPlateauDepolarisations.push_back(counter_of_plateau_depolarisations);
00174 
00175                     //update the counters.
00176                     ap_counter++;
00177                     ap_phase = BELOWTHRESHOLD;
00178 
00179                     //reinitialise counter of plateau depolarisations
00180                     counter_of_plateau_depolarisations = 0;
00181                 }
00182                 break;
00183         }
00184         prev_v = v;
00185         prev_t = t;
00186         prev_voltage_derivative = voltage_derivative;
00187     }
00188 
00189 
00190     // One last check. If the simulation ends halfway through an AP
00191     // i.e. if the vectors of onsets has more elements than the vectors
00192     // of peak and upstroke properties (that are updated at the end of the AP),
00193     // then we register the peak and upstroke values so far
00194     // for the last incomplete AP.
00195     if (mOnsets.size()>mMaxUpstrokeVelocities.size())
00196     {
00197         mMaxUpstrokeVelocities.push_back(max_upstroke_velocity);
00198         mPeakValues.push_back(current_peak);
00199         mTimesAtMaxUpstrokeVelocity.push_back(current_time_of_upstroke_velocity);
00200         mUnfinishedActionPotentials = true;
00201     }
00202 }
00203 
00204 
00205 std::vector<double> CellProperties::CalculateActionPotentialDurations(const double percentage)
00206 {
00207     CheckExceededThreshold();
00208 
00209     double prev_v = mrVoltage[0];
00210     unsigned APcounter=0;//will keep count of the APDs that we calculate
00211     bool apd_is_calculated=true;//this will ensure we hit the target only once per AP.
00212     std::vector<double> apds;
00213     double target = DBL_MAX;
00214     bool apd_starting_time_found=false;
00215     double apd_start_time=DBL_MAX;
00216 
00217     double t;
00218     double v;
00219     for (unsigned i=1; i<mrTime.size(); i++)
00220     {
00221         t = mrTime[i];
00222         v = mrVoltage[i];
00223 
00224         //First we make sure we stop calculating after the last AP has been calculated
00225         if (APcounter<mPeakValues.size())
00226         {
00227             //Set the target potential
00228             target = mRestingValues[APcounter]+0.01*(100-percentage)*(mPeakValues[APcounter]-mRestingValues[APcounter]);
00229 
00230             //if we reach the peak, we need to start to calculate an APD
00231             if (fabs(v-mPeakValues[APcounter])<=1e-6)
00232             {
00233                 apd_is_calculated = false;
00234             }
00235 
00236             // Start the timing where we first cross the target voltage
00237             if ( prev_v<v && prev_v<=target && v>=target && apd_starting_time_found==false)
00238             {
00239                 // Linear interpolation of target crossing time.
00240                 apd_start_time=t+( (target-prev_v)/(v-prev_v) )*(t-mrTime[i-1]);
00241                 apd_starting_time_found = true;
00242             }
00243 
00244             //if we hit the target while repolarising
00245             //and we are told this apd is not calculated yet.
00246             if ( prev_v>v && prev_v>=target && v<=target && apd_is_calculated==false)
00247             {
00248                 // Linear interpolation of target crossing time.
00249                 apds.push_back (t - apd_start_time + ( (target-prev_v)/(v-prev_v) )*(t-mrTime[i-1]) );
00250                 APcounter++;
00251                 apd_is_calculated = true;
00252                 apd_starting_time_found = false;
00253             }
00254         }
00255         prev_v = v;
00256     }
00257     if (apds.size() == 0)
00258     {
00259         EXCEPTION("No full action potential was recorded");
00260     }
00261     return apds;
00262 }
00263 
00264 std::vector<double> CellProperties::GetActionPotentialAmplitudes()
00265 {
00266     CheckExceededThreshold();
00267     unsigned size = mPeakValues.size();
00268     std::vector<double> amplitudes(size);
00269     for (unsigned i=0; i< size ;i++)
00270     {
00271         amplitudes[i] = (mPeakValues[i] - mRestingValues[i]);
00272     }
00273     return amplitudes;
00274 }
00275 
00276 void CellProperties::CheckExceededThreshold(void)
00277 {
00278     // mOnsets and mRestingValues are all
00279     // set at the same time, so checking one should suffice.
00280     if (mOnsets.empty())
00281     {
00282         // possible false error here if the simulation started at time < 0
00283         EXCEPTION("AP did not occur, never exceeded threshold voltage.");
00284     }
00285 }
00286 
00287 void CellProperties::CheckReturnedToThreshold(void)
00288 {
00289     // mPeakValues, mMaxUpstrokeVelocities and mTimesAtMaxUpstrokeVelocity are all
00290     // set at the same time so checking one should suffice.
00291     if (mMaxUpstrokeVelocities.empty())
00292     {
00293         EXCEPTION("AP did not occur, never descended past threshold voltage.");
00294     }
00295 }
00296 
00297 //
00298 // The Get std::vector<double> methods
00299 //
00300 
00301 std::vector<double> CellProperties::GetCycleLengths()
00302 {
00303     return mCycleLengths;
00304 }
00305 
00306 std::vector<double> CellProperties::GetRestingPotentials()
00307 {
00308     CheckExceededThreshold();
00309     return mRestingValues;
00310 }
00311 
00312 std::vector<double> CellProperties::GetPeakPotentials()
00313 {
00314     CheckReturnedToThreshold();
00315     return mPeakValues;
00316 }
00317 
00318 std::vector<double> CellProperties::GetMaxUpstrokeVelocities()
00319 {
00320     CheckReturnedToThreshold();
00321     return mMaxUpstrokeVelocities;
00322 }
00323 
00324 std::vector<double> CellProperties::GetTimesAtMaxUpstrokeVelocity()
00325 {
00326     CheckReturnedToThreshold();
00327     return mTimesAtMaxUpstrokeVelocity;
00328 }
00329 
00330 std::vector<double> CellProperties::GetAllActionPotentialDurations(const double percentage)
00331 {
00332     return CalculateActionPotentialDurations(percentage);
00333 }
00334 
00335 std::vector<unsigned> CellProperties::GetNumberOfAboveThresholdDepolarisationsForAllAps()
00336 {
00337     return mCounterOfPlateauDepolarisations;
00338 }
00339 //
00340 // The Get <double> methods
00341 //
00342 
00343 double CellProperties::GetLastPeakPotential()
00344 {
00345     CheckReturnedToThreshold();
00346     return mPeakValues.back();
00347 }
00348 
00349 double CellProperties::GetLastCompletePeakPotential()
00350 {
00351     CheckReturnedToThreshold();
00352     double peak_value;
00353     if (mUnfinishedActionPotentials && mPeakValues.size()>1u)
00354     {
00355         peak_value = mPeakValues[mPeakValues.size()-2u];
00356     }
00357     else if  (mUnfinishedActionPotentials && mPeakValues.size()==1u)
00358     {
00359         EXCEPTION("No peak potential matching a full action potential was recorded.");
00360     }
00361     else
00362     {
00363         peak_value = mPeakValues.back();
00364     }
00365     return peak_value;
00366 }
00367 
00368 double CellProperties::GetLastMaxUpstrokeVelocity()
00369 {
00370     CheckReturnedToThreshold();
00371     return mMaxUpstrokeVelocities.back();
00372 }
00373 
00374 double CellProperties::GetLastCompleteMaxUpstrokeVelocity()
00375 {
00376     CheckReturnedToThreshold();
00377     double max_upstroke;
00378     if (mUnfinishedActionPotentials && mMaxUpstrokeVelocities.size()>1u)
00379     {
00380         max_upstroke = mMaxUpstrokeVelocities[mMaxUpstrokeVelocities.size()-2u];
00381     }
00382     else if  (mUnfinishedActionPotentials && mMaxUpstrokeVelocities.size()==1u)
00383     {
00384         EXCEPTION("No MaxUpstrokeVelocity matching a full action potential was recorded.");
00385     }
00386     else
00387     {
00388         max_upstroke = mMaxUpstrokeVelocities.back();
00389     }
00390     return max_upstroke;
00391 }
00392 
00393 double CellProperties::GetTimeAtLastMaxUpstrokeVelocity()
00394 {
00395     CheckReturnedToThreshold();
00396     return mTimesAtMaxUpstrokeVelocity.back();
00397 }
00398 
00399 double CellProperties::GetTimeAtLastCompleteMaxUpstrokeVelocity()
00400 {
00401     CheckReturnedToThreshold();
00402     double max_upstroke_time;
00403     if (mUnfinishedActionPotentials && mTimesAtMaxUpstrokeVelocity.size()>1u)
00404     {
00405         max_upstroke_time = mTimesAtMaxUpstrokeVelocity[mTimesAtMaxUpstrokeVelocity.size()-2u];
00406     }
00407     else if  (mUnfinishedActionPotentials && mTimesAtMaxUpstrokeVelocity.size()==1u)
00408     {
00409         EXCEPTION("No TimeAtMaxUpstrokeVelocity matching a full action potential was recorded.");
00410     }
00411     else
00412     {
00413         max_upstroke_time = mTimesAtMaxUpstrokeVelocity.back();
00414     }
00415     return max_upstroke_time;
00416 }
00417 
00418 double CellProperties::GetLastActionPotentialDuration(const double percentage)
00419 {
00420     std::vector<double> apds = CalculateActionPotentialDurations(percentage);
00421     // We tested this returns a non-empty vector in the method.
00422     return apds.back();
00423 }
00424 
00425 unsigned CellProperties::GetNumberOfAboveThresholdDepolarisationsForLastAp()
00426 {
00427     return mCounterOfPlateauDepolarisations.back();
00428 }
00429 
00430 

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5