00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef _RUNANDCHECKIONICMODELS_HPP_
00031 #define _RUNANDCHECKIONICMODELS_HPP_
00032
00033 #include <vector>
00034 #include <string>
00035
00036 #include "OdeSolution.hpp"
00037
00038 #include "ColumnDataWriter.hpp"
00039 #include "ColumnDataReader.hpp"
00040
00041 #include "AbstractCardiacCell.hpp"
00042 #include "HeartConfig.hpp"
00043
00044 void RunOdeSolverWithIonicModel(AbstractCardiacCell* pOdeSystem,
00045 double endTime,
00046 std::string filename,
00047 int stepPerRow=100,
00048 bool doComputeExceptVoltage=true,
00049 bool useSamplingInterval=false);
00050
00051 void CheckCellModelResults(const std::string& rBaseResultsFilename,
00052 std::string validResultsBasename = "",
00053 double tolerance = 1e-3);
00054
00055 std::vector<double> GetVoltages(ColumnDataReader& rReader);
00056
00057 void CompareCellModelResults(std::string baseResultsFilename1, std::string baseResultsFilename2,
00058 double tolerance, bool vOnly=false, std::string folderName="TestIonicModels");
00059
00060
00061
00062
00063
00064
00065
00066
00067 #include <cmath>
00068
00069 void RunOdeSolverWithIonicModel(AbstractCardiacCell* pOdeSystem,
00070 double endTime,
00071 std::string filename,
00072 int stepPerRow,
00073 bool doComputeExceptVoltage,
00074 bool useSamplingInterval)
00075 {
00076 double start_time = 0.0;
00077
00078 if (doComputeExceptVoltage)
00079 {
00080
00081 std::vector<double> state_variables_copy = pOdeSystem->rGetStateVariables();
00082
00083
00084 double v_init = pOdeSystem->GetVoltage();
00085 pOdeSystem->ComputeExceptVoltage(start_time, start_time + 100.0 );
00086 double v_end = pOdeSystem->GetVoltage();
00087 TS_ASSERT_DELTA(v_init, v_end, 1e-6);
00088
00089
00090 pOdeSystem->SetVoltage(1e6);
00091 TS_ASSERT_DELTA(pOdeSystem->GetVoltage(), 1e6, 1e-6);
00092
00093
00094 pOdeSystem->SetStateVariables(state_variables_copy);
00095 }
00096
00097
00098 if (useSamplingInterval)
00099 {
00100 OdeSolution solution = pOdeSystem->Compute(start_time, endTime, HeartConfig::Instance()->GetOdeTimeStep() * stepPerRow);
00101 solution.WriteToFile("TestIonicModels", filename, "ms", 1, false, 4);
00102 }
00103 else
00104 {
00105 OdeSolution solution = pOdeSystem->Compute(start_time, endTime);
00106 solution.WriteToFile("TestIonicModels", filename, "ms", stepPerRow, false, 4);
00107 }
00108 }
00109
00110 std::vector<double> GetVoltages(ColumnDataReader& rReader)
00111 {
00112
00113
00114 std::vector<double> voltages;
00115 if (rReader.HasValues("V"))
00116 {
00117 voltages = rReader.GetValues("V");
00118 }
00119 else if (rReader.HasValues("membrane__V"))
00120 {
00121 voltages = rReader.GetValues("membrane__V");
00122 }
00123 else if (rReader.HasValues("membrane_voltage"))
00124 {
00125 voltages = rReader.GetValues("membrane_voltage");
00126 }
00127 else
00128 {
00129 EXCEPTION("Model membrane voltage not recognised.");
00130 }
00131 return voltages;
00132 }
00133
00134 std::vector<double> GetIntracellularCalcium(ColumnDataReader& rReader)
00135 {
00136
00137
00138 std::vector<double> cai;
00139 if (rReader.HasValues("CaI"))
00140 {
00141 cai = rReader.GetValues("CaI");
00142 }
00143 else if (rReader.HasValues("Cai"))
00144 {
00145 cai = rReader.GetValues("Cai");
00146 }
00147 else if (rReader.HasValues("cytosolic_calcium_concentration"))
00148 {
00149 cai = rReader.GetValues("cytosolic_calcium_concentration");
00150 }
00151 else
00152 {
00153 EXCEPTION("Model intracellular calcium is not recognised.");
00154 }
00155 return cai;
00156 }
00157
00158 std::vector<double> GetHGate(ColumnDataReader& rReader)
00159 {
00160
00161
00162 std::vector<double> h_values;
00163 if (rReader.HasValues("h"))
00164 {
00165 h_values = rReader.GetValues("h");
00166 }
00167 else if (rReader.HasValues("fast_sodium_current_h_gate__h"))
00168 {
00169 h_values = rReader.GetValues("fast_sodium_current_h_gate__h");
00170 }
00171 else
00172 {
00173 EXCEPTION("Model h gate is not recognised.");
00174 }
00175 return h_values;
00176 }
00177
00178
00179
00180
00181
00182 void CheckCellModelResults(const std::string& rBaseResultsFilename,
00183 std::string validResultsBasename,
00184 double tolerance)
00185 {
00186
00187 ColumnDataReader data_reader("TestIonicModels", rBaseResultsFilename);
00188 std::vector<double> times = data_reader.GetValues("Time");
00189 std::vector<double> voltages = GetVoltages(data_reader);
00190
00191 if (validResultsBasename == "")
00192 {
00193 validResultsBasename = rBaseResultsFilename;
00194 }
00195
00196 ColumnDataReader valid_reader("heart/test/data/ionicmodels", validResultsBasename + "ValidData",
00197 false);
00198 std::vector<double> valid_times = valid_reader.GetValues("Time");
00199 std::vector<double> valid_voltages = GetVoltages(valid_reader);
00200
00201 TS_ASSERT_EQUALS(times.size(), valid_times.size());
00202 for (unsigned i=0; i<valid_times.size(); i++)
00203 {
00204 TS_ASSERT_DELTA(times[i], valid_times[i], 1e-12);
00205
00206 TS_ASSERT_DELTA(voltages[i], valid_voltages[i], tolerance);
00207 }
00208 }
00209
00210 void CompareCellModelResults(std::string baseResultsFilename1, std::string baseResultsFilename2,
00211 double tolerance, bool vOnly, std::string folderName)
00212 {
00213
00214
00215 ColumnDataReader data_reader1(folderName, baseResultsFilename1);
00216 std::vector<double> times1 = data_reader1.GetValues("Time");
00217 std::vector<double> voltages1 = GetVoltages(data_reader1);
00218 std::vector<double> calcium1;
00219 std::vector<double> h1;
00220
00221 ColumnDataReader data_reader2(folderName, baseResultsFilename2);
00222 std::vector<double> times2 = data_reader2.GetValues("Time");
00223 std::vector<double> voltages2 = GetVoltages(data_reader2);
00224 std::vector<double> calcium2;
00225 std::vector<double> h2;
00226
00227 if (!vOnly)
00228 {
00229 calcium1 = GetIntracellularCalcium(data_reader1);
00230 h1 = GetHGate(data_reader1);
00231 calcium2 = GetIntracellularCalcium(data_reader2);
00232 h2 = GetHGate(data_reader2);
00233 }
00234
00235 TS_ASSERT(times1.size() >= times2.size());
00236 double last_v = voltages2[0];
00237 double tol = tolerance;
00238 for (unsigned i=0, j=0; i<times2.size(); i++)
00239 {
00240
00241 while (j<times1.size() && times1[j] < times2[i] - 1e-12)
00242 {
00243 j++;
00244 }
00245
00246
00247 if (fabs(voltages2[i] - last_v) > 0.05)
00248 {
00249 tol = tolerance * 25;
00250 }
00251 else
00252 {
00253 tol = tolerance;
00254 }
00255 last_v = voltages2[i];
00256
00257 TS_ASSERT_DELTA(times1[j], times2[i], 1e-12);
00258
00259 TS_ASSERT_DELTA(voltages1[j], voltages2[i], tol);
00260 if (!vOnly)
00261 {
00262 TS_ASSERT_DELTA(calcium1[j], calcium2[i], tol/100);
00263 TS_ASSERT_DELTA(h1[j], h2[i], tol/10);
00264 }
00265 }
00266 }
00267
00268
00269 #endif //_RUNANDCHECKIONICMODELS_HPP_