Chaste  Release::2018.1
NhsContractionModel.cpp
1 /*
2 
3 Copyright (c) 2005-2018, 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 #include "NhsContractionModel.hpp"
36 #include "OdeSystemInformation.hpp"
37 #include "EulerIvpOdeSolver.hpp"
38 #include "UblasCustomFunctions.hpp"
39 
40 #include <cmath>
41 
42 //
43 // Model-scope constant parameters
44 //
45 const double NhsContractionModel::mKon = 100;
46 const double NhsContractionModel::mKrefoff = 0.2;
47 const double NhsContractionModel::mGamma = 2;
49 const double NhsContractionModel::mAlphaR1 = 0.002;
50 const double NhsContractionModel::mAlphaR2 = 0.0017;
51 const double NhsContractionModel::mKZ = 0.15;
52 const unsigned NhsContractionModel::mNr = 3u;
53 const double NhsContractionModel::mBeta1 = -4;
54 const double NhsContractionModel::mAlpha0 = 0.008;
55 const unsigned NhsContractionModel::mN = 3u;
56 const double NhsContractionModel::mZp = 0.85;
57 const double NhsContractionModel::mCalcium50ref = 0.00105;
58 const double NhsContractionModel::mTref = 56.2;
59 const double NhsContractionModel::mBeta0 = 4.9;
60 const double NhsContractionModel::mA = 0.35;
61 const double NhsContractionModel::mA1 = -29;
62 const double NhsContractionModel::mA2 = 138;
63 const double NhsContractionModel::mA3 = 129;
64 const double NhsContractionModel::mAlpha1 = 0.03;
65 const double NhsContractionModel::mAlpha2 = 0.130;
66 const double NhsContractionModel::mAlpha3 = 0.625;
67 
68 
69 /*
70  * ============================== PRIVATE FUNCTIONS =====================================
71  */
73 {
74  double Ca50ref_times_one_plus_beta1_times_lam_minus_one = mCalcium50ref * (1 + mBeta1*(mLambda-1));
75  double one_plus_beta0_times_lam_minus_one_over_two_gamma = (1 + mBeta0*(mLambda-1))/(2*mGamma);
76 
77  mCalciumTrop50 = mCalciumTroponinMax * Ca50ref_times_one_plus_beta1_times_lam_minus_one;
78  mCalciumTrop50 /= (Ca50ref_times_one_plus_beta1_times_lam_minus_one + (1-one_plus_beta0_times_lam_minus_one_over_two_gamma)*mKrefoff/mKon);
79 }
80 
81 
83 {
84  double calcium_ratio_to_n = SmallPow(mCalciumTrop50/mCalciumTroponinMax, mN);
85 
86  double z_max = mAlpha0 - mK2*calcium_ratio_to_n;
87  z_max /= mAlpha0 + (mAlphaR1 + mK1)*calcium_ratio_to_n;
88 
89  return z * mTref * (1+mBeta0*(mLambda-1)) / z_max;
90 }
91 
92 /*
93  * ============================== PUBLIC FUNCTIONS =====================================
94  */
95 
97  : AbstractOdeBasedContractionModel(5) // five state variables
98 {
101 
102  mLambda = 1.0;
103  mDLambdaDt = 0.0;
104  mCalciumI = 0.0;
105 
106  // Initialise mCalciumTrop50!!
108 
109  double zp_to_n_plus_K_to_n = SmallPow(mZp,mNr) + SmallPow(mKZ,mNr);
110 
111  mK1 = mAlphaR2 * SmallPow(mZp,mNr-1) * mNr * SmallPow(mKZ,mNr);
112  mK1 /= zp_to_n_plus_K_to_n * zp_to_n_plus_K_to_n;
113 
114  mK2 = mAlphaR2 * SmallPow(mZp,mNr)/zp_to_n_plus_K_to_n;
115  mK2 *= 1 - mNr*SmallPow(mKZ,mNr)/zp_to_n_plus_K_to_n;
116 }
117 
118 void NhsContractionModel::SetStretchAndStretchRate(double lambda, double dlambdaDt)
119 {
120  assert(lambda>0.0);
121  mLambda = lambda;
122  mDLambdaDt = dlambdaDt;
123  // lambda changed so update mCalciumTrop50!!
125 }
126 
128 {
129  assert(rInputParameters.intracellularCalciumConcentration != DOUBLE_UNSET);
130  assert(rInputParameters.intracellularCalciumConcentration > 0.0);
131  mCalciumI = rInputParameters.intracellularCalciumConcentration;
132 }
133 
135 {
136  assert(calciumConcentration > 0.0);
137  mCalciumI = calciumConcentration;
138 }
139 
141 {
142  return mStateVariables[0];
143 }
144 
146  const std::vector<double> &rY,
147  std::vector<double> &rDY)
148 {
150 
151  const double& calcium_troponin = rY[0];
152  const double& z = rY[1];
153  const double& Q1 = rY[2];
154  const double& Q2 = rY[3];
155  const double& Q3 = rY[4];
156 
157  // check the state vars are in the expected range
158  // LCOV_EXCL_START
159  if (calcium_troponin < 0)
160  {
161  EXCEPTION("CalciumTrop concentration went negative");
162  }
163  if (z<0)
164  {
165  EXCEPTION("z went negative");
166  }
167  if (z>1)
168  {
169  EXCEPTION("z became greater than 1");
170  }
171  // LCOV_EXCL_STOP
172 
173 
174  double Q = Q1 + Q2 + Q3;
175  double T0 = CalculateT0(z);
176 
177  double Ta;
178  if (Q>0)
179  {
180  Ta = T0*(1+(2+mA)*Q)/(1+Q);
181  }
182  else
183  {
184  Ta = T0*(1+mA*Q)/(1-Q);
185  }
186 
187  rDY[0] = mKon * mCalciumI * ( mCalciumTroponinMax - calcium_troponin)
188  - mKrefoff * (1-Ta/(mGamma*mTref)) * calcium_troponin;
189 
190  double ca_trop_to_ca_trop50_ratio_to_n = SmallPow(calcium_troponin/mCalciumTrop50, mN);
191 
192  rDY[1] = mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n * (1-z)
193  - mAlphaR1 * z
194  - mAlphaR2 * SmallPow(z,mNr) / (SmallPow(z,mNr) + SmallPow(mKZ,mNr));
195 
196 
197  rDY[2] = mA1 * mDLambdaDt - mAlpha1 * Q1;
198  rDY[3] = mA2 * mDLambdaDt - mAlpha2 * Q2;
199  rDY[4] = mA3 * mDLambdaDt - mAlpha3 * Q3;
200 }
201 
202 
204 {
205  double T0 = CalculateT0(mStateVariables[1]);
206  double Q = mStateVariables[2]+mStateVariables[3]+mStateVariables[4];
207 
208  if (Q>0)
209  {
210  return T0*(1+(2+mA)*Q)/(1+Q);
211  }
212  else
213  {
214  return T0*(1+mA*Q)/(1-Q);
215  }
216 }
217 
218 template<>
220 {
221  this->mVariableNames.push_back("CalciumTroponin");
222  this->mVariableUnits.push_back("microMols");
223  this->mInitialConditions.push_back(0);
224 
225  this->mVariableNames.push_back("z");
226  this->mVariableUnits.push_back("");
227  this->mInitialConditions.push_back(0);
228 
229  this->mVariableNames.push_back("Q1");
230  this->mVariableUnits.push_back("");
231  this->mInitialConditions.push_back(0);
232 
233  this->mVariableNames.push_back("Q2");
234  this->mVariableUnits.push_back("");
235  this->mInitialConditions.push_back(0);
236 
237  this->mVariableNames.push_back("Q3");
238  this->mVariableUnits.push_back("");
239  this->mInitialConditions.push_back(0);
240 
241  this->mInitialised = true;
242 }
static boost::shared_ptr< OdeSystemInformation< ODE_SYSTEM > > Instance()
static const double mBeta0
void SetStretchAndStretchRate(double lambda, double dlambdaDt)
double CalculateT0(double z)
double SmallPow(double x, unsigned exponent)
static const double mTref
#define EXCEPTION(message)
Definition: Exception.hpp:143
static const double mKZ
static const unsigned mNr
static const double mZp
static const double mAlphaR1
const double DOUBLE_UNSET
Definition: Exception.hpp:56
static const double mAlpha0
void SetInputParameters(ContractionModelInputParameters &rInputParameters)
static const unsigned mN
static const double mAlpha3
static const double mCalciumTroponinMax
void SetIntracellularCalciumConcentration(double calciumConcentration)
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
static const double mBeta1
static const double mA
static const double mA2
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
static const double mCalcium50ref
static const double mA3
static const double mAlpha2
static const double mGamma
static const double mKon
static const double mAlpha1
static const double mAlphaR2
static const double mKrefoff
static const double mA1