PropagationPropertiesCalculator.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 "UblasIncludes.hpp"
00030 #include "PropagationPropertiesCalculator.hpp"
00031 #include "CellProperties.hpp"
00032 #include "Exception.hpp"
00033 #include <sstream>
00034 
00035 PropagationPropertiesCalculator::PropagationPropertiesCalculator(Hdf5DataReader* pDataReader,
00036         const std::string voltageName)
00037         : mpDataReader(pDataReader),
00038         mVoltageName(voltageName)
00039 {}
00040 
00041 PropagationPropertiesCalculator::~PropagationPropertiesCalculator()
00042 {
00043     // We don't own the data reader, so we don't destroy it.
00044 }
00045 
00046 double PropagationPropertiesCalculator::CalculateMaximumUpstrokeVelocity(unsigned globalNodeIndex)
00047 {
00048     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00049     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00050     CellProperties cell_props(voltages, times);
00051     return cell_props.GetLastMaxUpstrokeVelocity();
00052 }
00053 
00054 std::vector<double> PropagationPropertiesCalculator::CalculateAllMaximumUpstrokeVelocities(unsigned globalNodeIndex, double threshold)
00055 {
00056     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00057     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00058     CellProperties cell_props(voltages, times, threshold);
00059     return cell_props.GetMaxUpstrokeVelocities();
00060 }
00061 
00062 std::vector<double> PropagationPropertiesCalculator::CalculateUpstrokeTimes(unsigned globalNodeIndex, double threshold)
00063 {
00064     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00065     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00066     CellProperties cell_props(voltages, times, threshold);
00067     return cell_props.GetTimesAtMaxUpstrokeVelocity();
00068 }
00069 
00070 double PropagationPropertiesCalculator::CalculateActionPotentialDuration(const double percentage,
00071         unsigned globalNodeIndex)
00072 {
00073     if (percentage < 1.0 || percentage >= 100.0)
00074     {
00075         EXCEPTION("First argument of CalculateActionPotentialDuration() is expected to be a percentage");
00076     }
00077     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00078     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00079     CellProperties cell_props(voltages, times);
00080     return cell_props.GetLastActionPotentialDuration(percentage);
00081 }
00082 
00083 std::vector<double> PropagationPropertiesCalculator::CalculateAllActionPotentialDurations(const double percentage,
00084         unsigned globalNodeIndex, double threshold)
00085 {
00086     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00087     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00088     CellProperties cell_props(voltages, times, threshold);
00089     return cell_props.GetAllActionPotentialDurations(percentage);
00090 }
00091 
00092 double PropagationPropertiesCalculator::CalculatePeakMembranePotential(unsigned globalNodeIndex)
00093 {
00094     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00095     double max = -DBL_MAX;
00096     for(unsigned i=0; i<voltages.size(); i++)
00097     {
00098         if(voltages[i]>max)
00099         {
00100             max = voltages[i];
00101         }
00102     }
00103     return max;
00104 }
00105 
00106 double PropagationPropertiesCalculator::CalculateConductionVelocity(unsigned globalNearNodeIndex,
00107         unsigned globalFarNodeIndex,
00108         const double euclideanDistance)
00109 {
00110     double t_near = 0;
00111     double t_far = 0;
00112     std::vector<double> near_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNearNodeIndex);
00113     std::vector<double> far_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalFarNodeIndex);
00114     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00115 
00116     CellProperties near_cell_props(near_voltages, times);
00117     CellProperties far_cell_props(far_voltages, times);
00118 
00119     //The size of each vector is the number of APs that reached that node
00120     unsigned aps_near_node = near_cell_props.GetMaxUpstrokeVelocities().size();
00121     unsigned aps_far_node = far_cell_props.GetMaxUpstrokeVelocities().size();
00122 
00123     //These should never be empty. If so, an exception should have been thrown in the GetMaxUpstrokeVelocities() method.
00124     assert(aps_near_node > 0);
00125     assert(aps_far_node > 0);
00126 
00127     //if the same number of APs reached both nodes, get the last one...
00128     if (aps_near_node == aps_far_node)
00129     {
00130         t_near = near_cell_props.GetTimeAtLastMaxUpstrokeVelocity();
00131         t_far = far_cell_props.GetTimeAtLastMaxUpstrokeVelocity();
00132     }
00133     //...otherwise get the one with the smallest value, which is the last AP to reach both nodes
00134     //This prevents possible calculation of negative conduction velocities
00135     //for repeated stimuli
00136     else if (aps_near_node > aps_far_node)
00137     {
00138         t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_far_node-1];
00139         t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_far_node-1];
00140     }
00141     else
00142     {
00143         t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_near_node-1];
00144         t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_near_node-1];
00145     }
00146 
00147     if ((globalNearNodeIndex == globalFarNodeIndex) || ( fabs(t_far - t_near) < 1e-8))
00148     {
00149         // globalNearNodeIndex and globalFarNodeIndex are the same node, preventing a 0/0
00150         // or
00151         // AP number i is happening at the same time at nodes globalNearNodeIndex and globalFarNodeIndex
00152         return 0.0;
00153     }
00154     else
00155     {
00156         return euclideanDistance / (t_far - t_near);
00157     }
00158 
00159 
00160 }
00161 
00162 std::vector<double> PropagationPropertiesCalculator::CalculateAllConductionVelocities(unsigned globalNearNodeIndex,
00163                                                                 unsigned globalFarNodeIndex,
00164                                                                 const double euclideanDistance)
00165 {
00166     std::vector<double> conduction_velocities;
00167 
00168     std::vector<double> t_near;
00169     std::vector<double> t_far;
00170     unsigned number_of_aps = 0;
00171 
00172     std::vector<double> near_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNearNodeIndex);
00173     std::vector<double> far_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalFarNodeIndex);
00174     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00175 
00176     CellProperties near_cell_props(near_voltages, times);
00177     CellProperties far_cell_props(far_voltages, times);
00178 
00179     t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity();
00180     t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity();
00181 
00182     //exception should have been thrown within the GetTimesAtMaxUpstrokeVelocity method if the threshold is never reached
00183     //and these vectors are empty
00184     assert(t_near.size() !=0);
00185     assert(t_far.size() !=0);
00186 
00187     //Check the node where the least number of aps is reached.
00188     //We will calculate only where AP reached both nodes
00189     if (t_near.size() > t_far.size())
00190     {
00191         number_of_aps=t_far.size();
00192     }
00193     else
00194     {
00195        number_of_aps=t_near.size();
00196     }
00197     //now fill the vector
00198 
00199     if (globalNearNodeIndex == globalFarNodeIndex)
00200     {
00201         // globalNearNodeIndex and globalFarNodeIndex are the same node, preventing a 0/0
00202         for (unsigned i = 0 ; i < number_of_aps;i++)
00203         {
00204             conduction_velocities.push_back(0.0);
00205         }
00206     }
00207     else
00208     {
00209         for (unsigned i = 0 ; i < number_of_aps;i++)
00210         {
00211             if ( fabs(t_far[i] - t_near[i]) < 1e-8)
00212             {
00213                 // AP number i is happening at the same time at nodes globalNearNodeIndex and globalFarNodeIndex
00214                 conduction_velocities.push_back(0.0);
00215             }
00216             else
00217             {
00218                 conduction_velocities.push_back(euclideanDistance / (t_far[i] - t_near[i]));
00219             }
00220         }
00221     }
00222 
00223     return conduction_velocities;
00224 }
00225 
00226 std::vector<unsigned> PropagationPropertiesCalculator::CalculateAllAboveThresholdDepolarisations(unsigned globalNodeIndex, double threshold)
00227 {
00228     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00229     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00230     CellProperties cell_props(voltages, times, threshold);
00231     return cell_props.GetNumberOfAboveThresholdDepolarisationsForAllAps();
00232 }
00233 
00234 
00235 unsigned PropagationPropertiesCalculator::CalculateAboveThresholdDepolarisationsForLastAp(unsigned globalNodeIndex, double threshold)
00236 {
00237     std::vector<double> voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalNodeIndex);
00238     std::vector<double> times = mpDataReader->GetUnlimitedDimensionValues();
00239     CellProperties cell_props(voltages, times, threshold);
00240     return cell_props.GetNumberOfAboveThresholdDepolarisationsForLastAp();
00241 }

Generated by  doxygen 1.6.2