Chaste  Release::3.4
RunAndCheckIonicModels.hpp
1 /*
2 
3 Copyright (c) 2005-2016, 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 
37 #ifndef _RUNANDCHECKIONICMODELS_HPP_
38 #define _RUNANDCHECKIONICMODELS_HPP_
39 
40 #include <vector>
41 #include <string>
42 
43 #include "OdeSolution.hpp"
44 
45 #include "ColumnDataWriter.hpp"
46 #include "ColumnDataReader.hpp"
47 
48 #include "AbstractCardiacCell.hpp"
49 #include "HeartConfig.hpp"
50 
51 void RunOdeSolverWithIonicModel(AbstractCardiacCellInterface* pOdeSystem,
52  double endTime,
53  std::string filename,
54  int stepPerRow=100,
55  bool doComputeExceptVoltage=true,
56  bool useSamplingInterval=false);
57 
58 void CheckCellModelResults(const std::string& rBaseResultsFilename,
59  std::string validResultsBasename = "",
60  double tolerance = 1e-3,
61  std::string referenceFolder = "heart/test/data/ionicmodels");
62 
63 std::vector<double> GetVoltages(ColumnDataReader& rReader);
64 
65 void CompareCellModelResults(std::string baseResultsFilename1, std::string baseResultsFilename2,
66  double tolerance, bool vOnly=false, std::string folderName="TestIonicModels");
67 
68 
69 /*
70  * Note: we have to have the function definitions here, rather than in a .cpp file, if we're
71  * to include them in heart/src, as otherwise linking of the heart library fails because
72  * CxxTest routines are not defined.
73  */
74 
75 #include <cmath>
76 
77 void RunOdeSolverWithIonicModel(AbstractCardiacCellInterface* pOdeSystem,
78  double endTime,
79  std::string filename,
80  int stepPerRow,
81  bool doComputeExceptVoltage,
82  bool useSamplingInterval)
83 {
84  double start_time = 0.0;
85 
86  if (doComputeExceptVoltage)
87  {
88  // Store the current system state
89  std::vector<double> state_variables_copy = pOdeSystem->GetStdVecStateVariables();
90 
91  // Test ComputeExceptVoltage
92  double v_init = pOdeSystem->GetVoltage();
93  pOdeSystem->ComputeExceptVoltage(start_time, start_time + 100.0 /*ms*/);
94  double v_end = pOdeSystem->GetVoltage();
95  TS_ASSERT_DELTA(v_init, v_end, 1e-6);
96 
97  // Test SetVoltage
98  pOdeSystem->SetVoltage(1e6);
99  TS_ASSERT_DELTA(pOdeSystem->GetVoltage(), 1e6, 1e-6);
100 
101  // Reset the system
102  pOdeSystem->SetStateVariables(state_variables_copy);
103  }
104 
105  // Solve and write to file
106  if (useSamplingInterval)
107  {
108  OdeSolution solution = pOdeSystem->Compute(start_time, endTime, HeartConfig::Instance()->GetOdeTimeStep() * stepPerRow);
109  solution.WriteToFile("TestIonicModels", filename, "ms", 1, false, 4);
110  }
111  else
112  {
113  OdeSolution solution = pOdeSystem->Compute(start_time, endTime);
114  solution.WriteToFile("TestIonicModels", filename, "ms", stepPerRow, false, 4);
115  }
116 }
117 
118 std::vector<double> GetVoltages(ColumnDataReader& rReader)
119 {
120  //Rather Ugly, we can't currently guarantee what the name of the voltage column is,
121  //hence we try to cover the most common possibilities
122  std::vector<double> voltages;
123  if (rReader.HasValues("V"))
124  {
125  voltages = rReader.GetValues("V");
126  }
127  else if (rReader.HasValues("membrane_voltage"))
128  {
129  voltages = rReader.GetValues("membrane_voltage");
130  }
131  else if (rReader.HasValues("membrane__V"))
132  {
133  voltages = rReader.GetValues("membrane__V");
134  }
135  else
136  {
137  EXCEPTION("Model membrane voltage not recognised.");
138  }
139  return voltages;
140 }
141 
142 std::vector<double> GetIntracellularCalcium(ColumnDataReader& rReader)
143 {
144  //Rather Ugly, we can't currently guarantee what the name of the calcium column is,
145  //hence we try to cover the most common possibilities
146  std::vector<double> cai;
147  if (rReader.HasValues("CaI"))
148  {
149  cai = rReader.GetValues("CaI");
150  }
151  else if (rReader.HasValues("cytosolic_calcium_concentration"))
152  {
153  cai = rReader.GetValues("cytosolic_calcium_concentration");
154  }
155  else if (rReader.HasValues("Cai"))
156  {
157  cai = rReader.GetValues("Cai");
158  }
159  else
160  {
161  EXCEPTION("Model intracellular calcium is not recognised.");
162  }
163  return cai;
164 }
165 
166 std::vector<double> GetHGate(ColumnDataReader& rReader)
167 {
168  //Rather Ugly, we can't currently guarantee what the name of the h gate column is,
169  //hence we try to cover the most common possibilities
170  std::vector<double> h_values;
171  if (rReader.HasValues("h"))
172  {
173  h_values = rReader.GetValues("h");
174  }
175  else if (rReader.HasValues("fast_sodium_current_h_gate__h"))
176  {
177  h_values = rReader.GetValues("fast_sodium_current_h_gate__h");
178  }
179  else if (rReader.HasValues("membrane_fast_sodium_current_h_gate"))
180  {
181  h_values = rReader.GetValues("membrane_fast_sodium_current_h_gate");
182  }
183  else
184  {
185  EXCEPTION("Model h gate is not recognised.");
186  }
187  return h_values;
188 }
189 
190 /*
191  * Check the cell model against a previous version
192  * or another source e.g. Alan's COR
193  */
194 void CheckCellModelResults(const std::string& rBaseResultsFilename,
195  std::string validResultsBasename,
196  double tolerance,
197  std::string referenceFolder)
198 {
199  // read data entries for the new file
200  ColumnDataReader data_reader("TestIonicModels", rBaseResultsFilename);
201  std::vector<double> times = data_reader.GetValues("Time");
202  std::vector<double> voltages = GetVoltages(data_reader);
203 
204  if (validResultsBasename == "")
205  {
206  validResultsBasename = rBaseResultsFilename;
207  }
208 
209  ColumnDataReader valid_reader(referenceFolder, validResultsBasename + "ValidData",
210  false);
211  std::vector<double> valid_times = valid_reader.GetValues("Time");
212  std::vector<double> valid_voltages = GetVoltages(valid_reader);
213 
214  TS_ASSERT_EQUALS(times.size(), valid_times.size());
215  for (unsigned i=0; i<valid_times.size(); i++)
216  {
217  TS_ASSERT_DELTA(times[i], valid_times[i], 1e-12);
218  TS_ASSERT_DELTA(voltages[i], valid_voltages[i], tolerance);
219  }
220 }
221 
222 void CompareCellModelResults(std::string baseResultsFilename1, std::string baseResultsFilename2,
223  double tolerance, bool vOnly, std::string folderName)
224 {
225  // Compare 2 sets of results, e.g. from 2 different solvers for the same model.
226  // If the time series differ, the finer resolution must be given first.
227  ColumnDataReader data_reader1(folderName, baseResultsFilename1);
228  std::vector<double> times1 = data_reader1.GetValues("Time");
229  std::vector<double> voltages1 = GetVoltages(data_reader1);
230  std::vector<double> calcium1;
231  std::vector<double> h1;
232 
233  ColumnDataReader data_reader2(folderName, baseResultsFilename2);
234  std::vector<double> times2 = data_reader2.GetValues("Time");
235  std::vector<double> voltages2 = GetVoltages(data_reader2);
236  std::vector<double> calcium2;
237  std::vector<double> h2;
238 
239  if (!vOnly)
240  {
241  calcium1 = GetIntracellularCalcium(data_reader1);
242  h1 = GetHGate(data_reader1);
243  calcium2 = GetIntracellularCalcium(data_reader2);
244  h2 = GetHGate(data_reader2);
245  }
246 
247  TS_ASSERT(times1.size() >= times2.size());
248  double last_v = voltages2[0];
249  double tol = tolerance;
250  for (unsigned i=0, j=0; i<times2.size(); i++)
251  {
252  // Find corresponding time index
253  while (j<times1.size() && times1[j] < times2[i] - 1e-12)
254  {
255  j++;
256  }
257 
258  // Set tolerance higher in upstroke
259  if (fabs(voltages2[i] - last_v) > 0.05)
260  {
261  tol = tolerance * 25;
262  }
263  else
264  {
265  tol = tolerance;
266  }
267  last_v = voltages2[i];
268 
269  TS_ASSERT_DELTA(times1[j], times2[i], 1e-12);
270  // adjust tol to data
271  TS_ASSERT_DELTA(voltages1[j], voltages2[i], tol);
272  if (!vOnly)
273  {
274  TS_ASSERT_DELTA(calcium1[j], calcium2[i], tol/100);
275  TS_ASSERT_DELTA(h1[j], h2[i], tol/10);
276  }
277  }
278 }
279 
280 
281 #endif //_RUNANDCHECKIONICMODELS_HPP_
bool HasValues(const std::string &rVariableName)
virtual double GetVoltage()=0
#define EXCEPTION(message)
Definition: Exception.hpp:143
std::vector< double > GetValues(const std::string &rVariableName)
virtual std::vector< double > GetStdVecStateVariables()=0
virtual void SetStateVariables(const std::vector< double > &rVariables)=0
virtual OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)=0
void WriteToFile(std::string directoryName, std::string baseResultsFilename, std::string timeUnits, unsigned stepsPerRow=1, bool cleanDirectory=true, unsigned precision=8, bool includeDerivedQuantities=false)
virtual void ComputeExceptVoltage(double tStart, double tEnd)=0
virtual void SetVoltage(double voltage)=0
static HeartConfig * Instance()