Chaste  Release::3.4
TysonNovak2001OdeSystem.cpp
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 #include "TysonNovak2001OdeSystem.hpp"
37 #include "OdeSystemInformation.hpp"
38 #include "MathsCustomFunctions.hpp"
39 
40 TysonNovak2001OdeSystem::TysonNovak2001OdeSystem(std::vector<double> stateVariables)
42 {
44 
45  Init();
46 
47  if (stateVariables != std::vector<double>())
48  {
49  SetStateVariables(stateVariables);
50  }
51 }
52 
54 {
55  // Do nothing
56 }
57 
59 {
60  // Initialise model parameter values
61  mK1 = 0.04;
62  mK2d = 0.04;
63  mK2dd = 1.0;
64  mK2ddd = 1.0;
65  mCycB_threshold = 0.1;
66  mK3d = 1.0;
67  mK3dd = 10.0;
68  mK4d = 2.0;
69  mK4 = 35;
70  mJ3 = 0.04;
71  mJ4 = 0.04;
72  mK5d = 0.005;
73  mK5dd = 0.2;
74  mK6 = 0.1;
75  mJ5 = 0.3;
76  mN = 4u;
77  mK7 = 1.0;
78  mK8 = 0.5;
79  mJ7 = 1e-3;
80  mJ8 = 1e-3;
81  mMad = 1.0;
82  mK9 = 0.1;
83  mK10 = 0.02;
84  mK11 = 1.0;
85  mK12d = 0.2;
86  mK12dd = 50.0;
87  mK12ddd = 100.0;
88  mKeq = 1e3;
89  mK13 = 1.0;
90  mK14 = 1.0;
91  mK15d = 1.5;
92  mK15dd = 0.05;
93  mK16d = 1.0;
94  mK16dd = 3.0;
95  mJ15 = 0.01;
96  mJ16 = 0.01;
97  mMu = 0.01;
98  mMstar = 10.0;
99 }
100 
101 void TysonNovak2001OdeSystem::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
102 {
103  double x1 = rY[0];
104  double x2 = rY[1];
105  double x3 = rY[2];
106  double x4 = rY[3];
107  double x5 = rY[4];
108  double x6 = rY[5];
109 
110  double dx1 = 0.0;
111  double dx2 = 0.0;
112  double dx3 = 0.0;
113  double dx4 = 0.0;
114  double dx5 = 0.0;
115  double dx6 = 0.0;
116 
117  double temp1 = 0.0;
118  double temp2 = 0.0;
119  double temp3 = 0.0;
120 
130  dx1 = mK1-(mK2d+mK2dd*x2)*x1;
131 
132  // The commented line below models the start transition, no cycling, without Cdc20A
133 // temp1 = ((mK3d)*(1.0-x2))/(mJ3+1.0-x2);
134 
135  temp1 = ((mK3d+mK3dd*x4)*(1.0-x2))/(mJ3+1.0-x2);
136  temp2 = (mK4*x6*x1*x2)/(mJ4+x2);
137  dx2 = temp1-temp2;
138 
139  temp1 = mK5dd*(SmallPow(x1*x6/mJ5,mN)/(1+SmallPow(x1*x6/mJ5,mN)));
140  temp2 = mK6*x3;
141  dx3 = mK5d + temp1 - temp2;
142 
143  temp1 = (mK7*x5*(x3-x4))/(mJ7+x3-x4);
144  temp2 = (mK8*mMad*x4)/(mJ8+x4);
145  temp3 = mK6*x4;
146  dx4 = temp1 - temp2 - temp3;
147 
148  dx5 = mK9*x6*x1*(1.0-x5) - mK10*x5;
149 
150  dx6 = mMu*x6*(1.0-x6/mMstar);
151 
152  // Multiply by 60 beacuase the Tyson and Novak 2001 paper has time in minutes, not hours
153  rDY[0] = dx1*60.0;
154  rDY[1] = dx2*60.0;
155  rDY[2] = dx3*60.0;
156  rDY[3] = dx4*60.0;
157  rDY[4] = dx5*60.0;
158  rDY[5] = dx6*60.0;
159 }
160 
161 void TysonNovak2001OdeSystem::AnalyticJacobian(const std::vector<double>& rSolutionGuess, double** jacobian, double time, double timeStep)
162 {
163  timeStep *= 60.0; // to scale Jacobian so in hours not minutes
164  double x1 = rSolutionGuess[0];
165  double x2 = rSolutionGuess[1];
166  double x3 = rSolutionGuess[2];
167  double x4 = rSolutionGuess[3];
168  double x5 = rSolutionGuess[4];
169  double x6 = rSolutionGuess[5];
170 
171  // f1
172  double df1_dx1 = -mK2d - mK2dd*x2;
173  double df1_dx2 = -mK2dd*x1;
174 
175  jacobian[0][0] = 1-timeStep*df1_dx1;
176  jacobian[0][1] = -timeStep*df1_dx2;
177 
178  // f2
179  double df2_dx1 = -mK4*x6*x2/(mJ4+x2);
180  double df2_dx2 = -mJ3*(mK3d + mK3dd*x4)/(SmallPow((mJ3 + 1 - x2),2))
181  -mJ4*mK4*x6*x1/(SmallPow((mJ4+x2),2));
182  double df2_dx4 = mK3dd*(1-x2)/(mJ3+1-x2);
183  double df2_dx6 = -mK4*x1*x2/(mJ4+x2);
184 
185  jacobian[1][0] = -timeStep*df2_dx1;
186  jacobian[1][1] = 1-timeStep*df2_dx2;
187  jacobian[1][3] = -timeStep*df2_dx4;
188  jacobian[1][5] = -timeStep*df2_dx6;
189 
190  //f3
191  double z = x1*x6/mJ5;
192  double df3_dx1 = (mK5dd*x6/mJ5)*mN*SmallPow(z,mN-1)/(SmallPow((1-SmallPow(z,mN)),2));
193  double df3_dx3 = -mK6;
194  double df3_dx6 = (mK5dd*x1/mJ5)*mN*SmallPow(z,mN-1)/(SmallPow((1-SmallPow(z,mN)),2));
195 
196  jacobian[2][0] = -timeStep*df3_dx1;
197  jacobian[2][2] = 1-timeStep*df3_dx3;
198  jacobian[2][5] = -timeStep*df3_dx6;
199 
200  // f4
201  double df4_dx3 = mJ7*mK7*x5/(SmallPow(mJ7+x3-x4,2));
202  double df4_dx4 = -mJ7*mK7*x5/(SmallPow(mJ7+x3-x4,2)) - mK6 - mJ8*mK8*mMad/(SmallPow(mJ8+x4,2));
203  double df4_dx5 = mK7*(x3-x4)/(mJ7+x3-x4);
204 
205  jacobian[3][2] = -timeStep*df4_dx3;
206  jacobian[3][3] = 1-timeStep*df4_dx4;
207  jacobian[3][4] = -timeStep*df4_dx5;
208 
209  // f5
210  double df5_dx1 = mK9*x6*(1-x5);
211  double df5_dx5 = -mK10 - mK9*x6*x1;
212  double df5_dx6 = mK9*x1*(1-x5);
213 
214  jacobian[4][0] = -timeStep*df5_dx1;
215  jacobian[4][4] = 1-timeStep*df5_dx5;
216  jacobian[4][5] = -timeStep*df5_dx6;
217 
218  // f6
219  double df6_dx6 = mMu - 2*mMu*x6/mMstar;
220 
221  jacobian[5][5] = 1-timeStep*df6_dx6;
222 }
223 
224 bool TysonNovak2001OdeSystem::CalculateStoppingEvent(double time, const std::vector<double>& rY)
225 {
226  std::vector<double> dy(rY.size());
227  EvaluateYDerivatives(time, rY, dy);
228 
229  // Only call this a stopping condition if the mass of the cell is over 0.6
230  // (normally cycles from 0.5-1.0 ish!)
231  return ( (rY[5] > 0.6 ) && (rY[0] < mCycB_threshold) && dy[0] < 0.0 );
232 }
233 
234 double TysonNovak2001OdeSystem::CalculateRootFunction(double time, const std::vector<double>& rY)
235 {
236  std::vector<double> dy(rY.size());
237  EvaluateYDerivatives(time, rY, dy);
238 
239  // Only call this a stopping condition if the mass of the cell is over 0.6
240  // (normally cycles from 0.5-1.0 ish!)
241  if (rY[5]<0.6)
242  {
243  return 1.0;
244  }
245 
246  if (dy[0] >= 0.0)
247  {
248  return 1.0;
249  }
250  return rY[0]-mCycB_threshold;
251 }
252 
253 template<>
255 {
256  /*
257  * Initialise state variables.
258  *
259  * These initial conditions are the approximate steady state
260  * solution values while the commented out conditions are taken
261  * from the Tyson and Novak 2001 paper.
262  */
263  this->mVariableNames.push_back("CycB");
264  this->mVariableUnits.push_back("nM");
265 // this->mInitialConditions.push_back(0.1);
266  this->mInitialConditions.push_back(0.099999999999977);
267 
268  this->mVariableNames.push_back("Cdh1");
269  this->mVariableUnits.push_back("nM");
270 // this->mInitialConditions.push_back(9.8770e-01);
271  this->mInitialConditions.push_back(0.989026454281841);
272 
273  this->mVariableNames.push_back("Cdc20T");
274  this->mVariableUnits.push_back("nM");
275 // this->mInitialConditions.push_back(1.5011e+00);
276  this->mInitialConditions.push_back(1.547942029285891);
277 
278  this->mVariableNames.push_back("Cdc20A");
279  this->mVariableUnits.push_back("nM");
280 // this->mInitialConditions.push_back(1.2924e+00);
281  this->mInitialConditions.push_back(1.421110920155839);
282 
283  this->mVariableNames.push_back("IEP");
284  this->mVariableUnits.push_back("nM");
285 // this->mInitialConditions.push_back(6.5405e-01);
286  this->mInitialConditions.push_back(0.672838844290094);
287 
288  this->mVariableNames.push_back("mass");
289  this->mVariableUnits.push_back("");
290 // this->mInitialConditions.push_back(4.7039e-01);
291  this->mInitialConditions.push_back(0.970831277863956 / 2);
292 
293  this->mInitialised = true;
294 }
295 
296 // Serialization for Boost >= 1.36
static boost::shared_ptr< OdeSystemInformation< ODE_SYSTEM > > Instance()
double SmallPow(double x, unsigned exponent)
virtual void AnalyticJacobian(const std::vector< double > &rSolutionGuess, double **jacobian, double time, double timeStep)
TysonNovak2001OdeSystem(std::vector< double > stateVariables=std::vector< double >())
void SetStateVariables(const std::vector< double > &rStateVariables)
double CalculateRootFunction(double time, const std::vector< double > &rY)
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
bool CalculateStoppingEvent(double time, const std::vector< double > &rY)
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
#define CHASTE_CLASS_EXPORT(T)