Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
PropagationPropertiesCalculator.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 "UblasIncludes.hpp"
37#include "PropagationPropertiesCalculator.hpp"
38#include "CellProperties.hpp"
39#include "Exception.hpp"
40#include <sstream>
41#include "HeartEventHandler.hpp"
42
44 const std::string voltageName)
45 : mpDataReader(pDataReader),
46 mVoltageName(voltageName),
47 mTimes(mpDataReader->GetUnlimitedDimensionValues()),
48 mCachedNodeGlobalIndex(UNSIGNED_UNSET)
49{}
50
52{
53 // We don't own the data reader, so we don't destroy it.
54}
55
57{
58 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
59 CellProperties cell_props(r_voltages, mTimes);
60 return cell_props.GetLastMaxUpstrokeVelocity();
61}
62
63std::vector<double> PropagationPropertiesCalculator::CalculateAllMaximumUpstrokeVelocities(unsigned globalNodeIndex, double threshold)
64{
65 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
66 CellProperties cell_props(r_voltages, mTimes, threshold);
67 return cell_props.GetMaxUpstrokeVelocities();
68}
69
70std::vector<double> PropagationPropertiesCalculator::CalculateUpstrokeTimes(unsigned globalNodeIndex, double threshold)
71{
72 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
73 CellProperties cell_props(r_voltages, mTimes, threshold);
74 return cell_props.GetTimesAtMaxUpstrokeVelocity();
75}
76
78 unsigned globalNodeIndex)
79{
80 if (percentage < 1.0 || percentage >= 100.0)
81 {
82 EXCEPTION("First argument of CalculateActionPotentialDuration() is expected to be a percentage");
83 }
84 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
85 CellProperties cell_props(r_voltages, mTimes);
86 return cell_props.GetLastActionPotentialDuration(percentage);
87}
88
90 unsigned globalNodeIndex,
91 double threshold)
92{
93 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
94 CellProperties cell_props(r_voltages, mTimes, threshold);
95 return cell_props.GetAllActionPotentialDurations(percentage);
96}
97
99 const double percentage,
100 unsigned lowerNodeIndex,
101 unsigned upperNodeIndex,
102 double threshold)
103{
104 std::vector<std::vector<double> > output_data;
105 output_data.reserve(upperNodeIndex-lowerNodeIndex+1);
106 unsigned num_nodes_per_data_block = 100; // number of nodes
107 unsigned num_complete_blocks = (upperNodeIndex-lowerNodeIndex) / num_nodes_per_data_block;
108 unsigned size_last_block = (upperNodeIndex-lowerNodeIndex) % num_nodes_per_data_block;
109
110 for (unsigned block_num=0;
111 block_num<num_complete_blocks+1;
112 block_num++)
113 {
114 unsigned num_nodes_to_read;
115 if (block_num != num_complete_blocks)
116 {
117 num_nodes_to_read = num_nodes_per_data_block;
118 }
119 else
120 {
121 num_nodes_to_read = size_last_block;
122 }
123
124 if (num_nodes_to_read > 0)
125 {
126 // Read a big block of data
127 unsigned low_node = lowerNodeIndex + block_num*num_nodes_per_data_block;
128 unsigned high_node = low_node + num_nodes_to_read;
129 std::vector<std::vector<double> > voltages = mpDataReader->GetVariableOverTimeOverMultipleNodes(mVoltageName, low_node, high_node);
130
131 for (unsigned node_within_block=0;
132 node_within_block < num_nodes_to_read;
133 node_within_block++)
134 {
135 std::vector<double>& r_voltages = voltages[node_within_block];
136 CellProperties cell_props(r_voltages, mTimes, threshold);
137 std::vector<double> apds;
138 try
139 {
140 apds = cell_props.GetAllActionPotentialDurations(percentage);
141 assert(apds.size() != 0);
142 }
143 catch (Exception& e)
144 {
145 assert(e.GetShortMessage()=="No full action potential was recorded" ||
146 e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage.");
147 apds.push_back(0);
148 assert(apds.size() == 1);
149 }
150 output_data.push_back(apds);
151 }
152 }
153 }
154 return output_data;
155}
156
158{
159 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
160 double max = -DBL_MAX;
161 for (unsigned i=0; i<r_voltages.size(); i++)
162 {
163 if (r_voltages[i]>max)
164 {
165 max = r_voltages[i];
166 }
167 }
168 return max;
169}
170
172 unsigned globalFarNodeIndex,
173 const double euclideanDistance)
174{
175 double t_near = 0;
176 double t_far = 0;
177 std::vector<double>& r_near_voltages = rGetCachedVoltages(globalNearNodeIndex);
178 std::vector<double> far_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalFarNodeIndex);
179
180 CellProperties near_cell_props(r_near_voltages, mTimes);
181 CellProperties far_cell_props(far_voltages, mTimes);
182
183 //The size of each vector is the number of APs that reached that node
184 unsigned aps_near_node = near_cell_props.GetMaxUpstrokeVelocities().size();
185 unsigned aps_far_node = far_cell_props.GetMaxUpstrokeVelocities().size();
186
187 //These should never be empty. If so, an exception should have been thrown in the GetMaxUpstrokeVelocities() method.
188 assert(aps_near_node > 0);
189 assert(aps_far_node > 0);
190
191 //if the same number of APs reached both nodes, get the last one...
192 if (aps_near_node == aps_far_node)
193 {
194 t_near = near_cell_props.GetTimeAtLastMaxUpstrokeVelocity();
195 t_far = far_cell_props.GetTimeAtLastMaxUpstrokeVelocity();
196 }
197 //...otherwise get the one with the smallest value, which is the last AP to reach both nodes
198 //This prevents possible calculation of negative conduction velocities
199 //for repeated stimuli
200 else if (aps_near_node > aps_far_node)
201 {
202 t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_far_node-1];
203 t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_far_node-1];
204 }
205 else
206 {
207 t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_near_node-1];
208 t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity()[aps_near_node-1];
209 }
210
212 if ((globalNearNodeIndex == globalFarNodeIndex) || ( fabs(t_far - t_near) < 1e-8))
213 {
214 // globalNearNodeIndex and globalFarNodeIndex are the same node, preventing a 0/0
215 // or
216 // AP number i is happening at the same time at nodes globalNearNodeIndex and globalFarNodeIndex
217 return 0.0;
218 }
219 else
220 {
221 return euclideanDistance / (t_far - t_near);
222 }
223
224
225}
226
227std::vector<double> PropagationPropertiesCalculator::CalculateAllConductionVelocities(unsigned globalNearNodeIndex,
228 unsigned globalFarNodeIndex,
229 const double euclideanDistance)
230{
231 std::vector<double> conduction_velocities;
232
233 std::vector<double> t_near;
234 std::vector<double> t_far;
235 unsigned number_of_aps = 0;
236
237 std::vector<double>& r_near_voltages = rGetCachedVoltages(globalNearNodeIndex);
238 std::vector<double> far_voltages = mpDataReader->GetVariableOverTime(mVoltageName, globalFarNodeIndex);
239
240 CellProperties near_cell_props(r_near_voltages, mTimes);
241 CellProperties far_cell_props(far_voltages, mTimes);
242
243 t_near = near_cell_props.GetTimesAtMaxUpstrokeVelocity();
244 t_far = far_cell_props.GetTimesAtMaxUpstrokeVelocity();
245
246 //exception should have been thrown within the GetTimesAtMaxUpstrokeVelocity method if the threshold is never reached
247 //and these vectors are empty
248 assert(t_near.size() !=0);
249 assert(t_far.size() !=0);
250
251 //Check the node where the least number of aps is reached.
252 //We will calculate only where AP reached both nodes
253 if (t_near.size() > t_far.size())
254 {
255 number_of_aps = t_far.size();
256 }
257 else
258 {
259 number_of_aps = t_near.size();
260 }
261 //now fill the vector
262
263 if (globalNearNodeIndex == globalFarNodeIndex)
264 {
265 // globalNearNodeIndex and globalFarNodeIndex are the same node, preventing a 0/0
266 for (unsigned i = 0 ; i < number_of_aps;i++)
267 {
268 conduction_velocities.push_back(0.0);
269 }
270 }
271 else
272 {
273 for (unsigned i = 0 ; i < number_of_aps;i++)
274 {
276 if (fabs(t_far[i] - t_near[i]) < 1e-8)
277 {
278 // AP number i is happening at the same time at nodes globalNearNodeIndex and globalFarNodeIndex
279 conduction_velocities.push_back(0.0);
280 }
281 else
282 {
283 conduction_velocities.push_back(euclideanDistance / (t_far[i] - t_near[i]));
284 }
285 }
286 }
287
288 return conduction_velocities;
289}
290
291
292std::vector<unsigned> PropagationPropertiesCalculator::CalculateAllAboveThresholdDepolarisations(unsigned globalNodeIndex, double threshold)
293{
294 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
295 CellProperties cell_props(r_voltages, mTimes, threshold);
297}
298
299
301{
302 std::vector<double>& r_voltages = rGetCachedVoltages(globalNodeIndex);
303 CellProperties cell_props(r_voltages, mTimes, threshold);
305}
306
307
308std::vector<double>& PropagationPropertiesCalculator::rGetCachedVoltages(unsigned globalNodeIndex)
309{
310 if (globalNodeIndex != mCachedNodeGlobalIndex)
311 {
313 mCachedNodeGlobalIndex = globalNodeIndex;
314 }
315 return mCachedVoltages;
316}
317
322
323
#define EXCEPTION(message)
const unsigned UNSIGNED_UNSET
Definition Exception.hpp:53
std::vector< double > GetTimesAtMaxUpstrokeVelocity()
unsigned GetNumberOfAboveThresholdDepolarisationsForLastAp()
std::vector< unsigned > GetNumberOfAboveThresholdDepolarisationsForAllAps()
std::vector< double > GetAllActionPotentialDurations(const double percentage)
std::vector< double > GetMaxUpstrokeVelocities()
double GetTimeAtLastMaxUpstrokeVelocity()
double GetLastMaxUpstrokeVelocity()
double GetLastActionPotentialDuration(const double percentage)
std::vector< std::vector< double > > GetVariableOverTimeOverMultipleNodes(const std::string &rVariableName, unsigned lowerIndex, unsigned upperIndex)
std::vector< double > GetVariableOverTime(const std::string &rVariableName, unsigned nodeIndex)
double CalculateConductionVelocity(unsigned globalNearNodeIndex, unsigned globalFarNodeIndex, const double euclideanDistance)
double CalculateMaximumUpstrokeVelocity(unsigned globalNodeIndex)
double CalculatePeakMembranePotential(unsigned globalNodeIndex)
PropagationPropertiesCalculator(Hdf5DataReader *pDataReader, const std::string voltageName="V")
unsigned CalculateAboveThresholdDepolarisationsForLastAp(unsigned globalNodeIndex, double threshold)
std::vector< unsigned > CalculateAllAboveThresholdDepolarisations(unsigned globalNodeIndex, double threshold)
void SetHdf5DataReader(Hdf5DataReader *pDataReader)
std::vector< std::vector< double > > CalculateAllActionPotentialDurationsForNodeRange(const double percentage, unsigned lowerNodeIndex, unsigned upperNodeIndex, double threshold)
std::vector< double > CalculateAllConductionVelocities(unsigned globalNearNodeIndex, unsigned globalFarNodeIndex, const double euclideanDistance)
double CalculateActionPotentialDuration(const double percentage, unsigned globalNodeIndex)
std::vector< double > CalculateUpstrokeTimes(unsigned globalNodeIndex, double threshold)
std::vector< double > CalculateAllActionPotentialDurations(const double percentage, unsigned globalNodeIndex, double threshold)
std::vector< double > CalculateAllMaximumUpstrokeVelocities(unsigned globalNodeIndex, double threshold)
std::vector< double > & rGetCachedVoltages(unsigned globalNodeIndex)