Chaste  Release::3.4
WntCellCycleOdeSystem.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 "WntCellCycleOdeSystem.hpp"
37 #include "CellwiseOdeSystemInformation.hpp"
38 #include "IsNan.hpp"
39 
40 // These #includes are needed for the constructor and EvaluateYDerivatives()
41 #include "ApcOneHitCellMutationState.hpp"
42 #include "ApcTwoHitCellMutationState.hpp"
43 #include "BetaCateninOneHitCellMutationState.hpp"
44 
46  boost::shared_ptr<AbstractCellMutationState> pMutationState,
47  std::vector<double> stateVariables)
48  : AbstractOdeSystem(9),
49  mpMutationState(pMutationState),
50  mWntLevel(wntLevel)
51 {
53 
68  Init(); // set up parameter values
69 
70  // Set up a Wnt signalling pathway in a steady state
71  double destruction_level = ma5d/(ma4d*wntLevel+ma5d);
72  double beta_cat_level_1 = -1.0;
73  double beta_cat_level_2 = -1.0;
74 
75  if (!mpMutationState)
76  {
77  // No mutations specified
78  }
80  {
81  // APC +/- : only half are active
82  beta_cat_level_1 = 0.5*ma2d/(ma2d+0.5*ma3d*destruction_level);
83  beta_cat_level_2 = 0.5*ma2d/(ma2d+0.5*ma3d*destruction_level);
84  }
86  {
87  // APC -/-
88  destruction_level = 0.0; // no active destruction complex
89  beta_cat_level_1 = 0.5; // fully active beta-catenin
90  beta_cat_level_2 = 0.5; // fully active beta-catenin
91  }
93  {
94  // Beta-cat delta 45
95  beta_cat_level_1 = 0.5*ma2d/(ma2d+ma3d*destruction_level);
96  beta_cat_level_2 = 0.5;
97  }
98  else
99  {
100  // healthy cells
101  beta_cat_level_1 = 0.5*ma2d/(ma2d+ma3d*destruction_level);
102  beta_cat_level_2 = 0.5*ma2d/(ma2d+ma3d*destruction_level);
103  }
104 
105  // Cell-specific initial conditions
106  SetDefaultInitialCondition(5, destruction_level);
107  SetDefaultInitialCondition(6, beta_cat_level_1);
108  SetDefaultInitialCondition(7, beta_cat_level_2);
109  SetDefaultInitialCondition(8, wntLevel);
110 
111  if (stateVariables != std::vector<double>())
112  {
113  SetStateVariables(stateVariables);
114  }
115 }
116 
117 void WntCellCycleOdeSystem::SetMutationState(boost::shared_ptr<AbstractCellMutationState> pMutationState)
118 {
119  mpMutationState = pMutationState;
120 }
121 
123 {
124  // Do nothing
125 }
126 
128 {
129  // Initialise model parameter values
130  // Swat (2004) Parameters
131  double k1 = 1.0;
132  double k2 = 1.6;
133  double k3 = 0.05;
134  double k16 = 0.4;
135  double k34 = 0.04;
136  double k43 = 0.01;
137  double k61 = 0.3;
138  double k23 = 0.3;
139  double a = 0.04;
140  double J11 = 0.5;
141  double J12 = 5.0;
142  double J61 = 5.0;
143  double J62 = 8.0;
144  double J13 = 0.002;
145  double J63 = 2.0;
146  double Km1 = 0.5;
147  double Km2 = 4.0;
148  double Km4 = 0.3;
149  double kp = 0.05;
150  double phi_pRb = 0.005;
151  double phi_E2F1 = 0.1;
152  double phi_CycDi = 0.023;
153  double phi_CycDa = 0.03;
154  double phi_pRbp = 0.06;
155 
156  // Mirams et al. parameter values
157  double a1 = 0.423;
158  double a2 = 2.57e-4;
159  double a3 = 1.72;
160  double a4 = 10.0;
161  double a5 = 0.5;
162  double WntMax = 10.0;
163  double mitogenic_factorF = 6.0e-4;
164  double APC_Total = 0.02;
165 
166  // Non-dimensionalise...
167  mk2d = k2/(Km2*phi_E2F1);
168  mk3d = k3*a1*mitogenic_factorF/(Km4*phi_E2F1*a2);
169  mk34d = k34/phi_E2F1;
170  mk43d = k43/phi_E2F1;
171  mk23d = k23*Km2/(Km4*phi_E2F1);
172  mad = a/Km2;
173  mJ11d = J11*phi_E2F1/k1;
174  mJ12d = J12*phi_E2F1/k1;
175  mJ13d = J13*phi_E2F1/k1;
176  mJ61d = J61*phi_E2F1/k1;
177  mJ62d = J62*phi_E2F1/k1;
178  mJ63d = J63*phi_E2F1/k1;
179  mKm1d = Km1/Km2;
180  mkpd = kp/(Km2*phi_E2F1);
181  mphi_r = phi_pRb/phi_E2F1;
182  mphi_i = phi_CycDi/phi_E2F1;
183  mphi_j = phi_CycDa/phi_E2F1;
184  mphi_p = phi_pRbp/phi_E2F1;
185  ma2d = a2/phi_E2F1;
186  ma3d = a3*APC_Total/phi_E2F1;
187  ma4d = a4*WntMax/phi_E2F1;
188  ma5d = a5/phi_E2F1;
189  mk16d = k16*Km4/phi_E2F1;
190  mk61d = k61/phi_E2F1;
191  mPhiE2F1 = phi_E2F1;
192 }
193 
194 void WntCellCycleOdeSystem::EvaluateYDerivatives(double time, const std::vector<double>& rY, std::vector<double>& rDY)
195 {
196  double r = rY[0];
197  double e = rY[1];
198  double i = rY[2];
199  double j = rY[3];
200  double p = rY[4];
201  double c = rY[5];
202  double b1 = rY[6];
203  double b2 = rY[7];
204  double wnt_level = rY[8];
205 
206  double dx1 = 0.0;
207  double dx2 = 0.0;
208  double dx3 = 0.0;
209  double dx4 = 0.0;
210  double dx5 = 0.0;
211  double dx6 = 0.0;
212  double dx7 = 0.0;
213  double dx8 = 0.0;
214 
215  /*
216  * The variables are
217  * 1. r = pRb
218  * 2. e = E2F1
219  * 3. i = CycD (inactive)
220  * 4. j = CycD (active)
221  * 5. p = pRb-p
222  * 6. c = APC (Active)
223  * 7. b = Beta-Catenin
224  */
225 
226  // Bit back-to-front, but work out the Wnt section first...
227 
228  // Mutations take effect by altering the level of beta-catenin
229  if (!mpMutationState)
230  {
231  // No mutations specified
232  }
233  else if (mpMutationState->IsType<ApcOneHitCellMutationState>()) // APC +/-
234  {
235  dx6 = ma5d*(1.0-c) - ma4d*wnt_level*c;
236  dx7 = ma2d*(0.5-b1) - 0.5*ma3d*b1*c;
237  dx8 = ma2d*(0.5-b2) - 0.5*ma3d*b2*c;
238  }
239  else if (mpMutationState->IsType<ApcTwoHitCellMutationState>()) // APC -/-
240  {
241  dx6 = 0.0;
242  dx7 = ma2d*(0.5-b1);
243  dx8 = ma2d*(0.5-b2);
244  }
245  else if (mpMutationState->IsType<BetaCateninOneHitCellMutationState>()) // Beta-Cat D45
246  {
247  dx6 = ma5d*(1.0-c) - ma4d*wnt_level*c;
248  dx7 = ma2d*(0.5-b1) - ma3d*b1*c;
249  dx8 = ma2d*(0.5-b2);
250  }
251  else
252  {
253  // da
254  dx6 = ma5d*(1.0-c) - ma4d*wnt_level*c;
255  // db
256  dx7 = ma2d*(0.5-b1) - ma3d*b1*c;
257  dx8 = ma2d*(0.5-b2) - ma3d*b2*c;
258  }
259 
260  // Now the cell cycle stuff...
261 
262  // dr
263  dx1 = e/(mKm1d+e)*mJ11d/(mJ11d+r)*mJ61d/(mJ61d+p) - mk16d*r*j+mk61d*p-mphi_r*r;
264  // de
265  dx2 = mkpd+mk2d*(mad*mad+e*e)/(1+e*e)*mJ12d/(mJ12d+r)*mJ62d/(mJ62d+p) - e;
266  // di
267  dx3 = mk3d*(b1+b2) + mk23d*e*mJ13d/(mJ13d+r)*mJ63d/(mJ63d+p) + mk43d*j - mk34d*i*j/(1+j) - mphi_i*i;
268  // dj
269  dx4 = mk34d*i*j/(1+j) - (mk43d+mphi_j)*j;
270  // dp
271  dx5 = mk16d*r*j - mk61d*p - mphi_p*p;
272 
273  double factor = mPhiE2F1*60.0; // convert non-dimensional d/dt s to d/dt in hours
274 
275  rDY[0] = dx1*factor;
276  rDY[1] = dx2*factor;
277  rDY[2] = dx3*factor;
278  rDY[3] = dx4*factor;
279  rDY[4] = dx5*factor;
280  rDY[5] = dx6*factor;
281  rDY[6] = dx7*factor; // beta-cat allele 1
282  rDY[7] = dx8*factor; // beta-cat allele 2
283  rDY[8] = 0.0; // do not change the Wnt level
284 }
285 
286 const boost::shared_ptr<AbstractCellMutationState> WntCellCycleOdeSystem::GetMutationState() const
287 {
288  return mpMutationState;
289 }
290 
291 bool WntCellCycleOdeSystem::CalculateStoppingEvent(double time, const std::vector<double>& rY)
292 {
293  double r = rY[0];
294  double e = rY[1];
295  double p = rY[4];
296  double dY1 = mkpd+mk2d*(mad*mad+e*e)/(1+e*e)*mJ12d/(mJ12d+r)*mJ62d/(mJ62d+p) - e;
297  double factor = mPhiE2F1*60.0; // Convert non-dimensional d/dt s to d/dt in hours.
298  dY1 = dY1*factor;
299 
300  assert(!std::isnan(rY[1]));
301  assert(!std::isnan(dY1));
302  return (rY[1] > 1.0 && dY1 > 0.0);
303 }
304 
305 double WntCellCycleOdeSystem::CalculateRootFunction(double time, const std::vector<double>& rY)
306 {
307  return rY[1] - 1.0;
308 }
309 
310 template<>
312 {
313  this->mVariableNames.push_back("pRb");
314  this->mVariableUnits.push_back("non_dim");
315  this->mInitialConditions.push_back(7.357000000000000e-01);
316 
317  this->mVariableNames.push_back("E2F1");
318  this->mVariableUnits.push_back("non_dim");
319  this->mInitialConditions.push_back(1.713000000000000e-01);
320 
321  this->mVariableNames.push_back("CycD_i");
322  this->mVariableUnits.push_back("non_dim");
323  this->mInitialConditions.push_back(6.900000000000001e-02);
324 
325  this->mVariableNames.push_back("CycD_a");
326  this->mVariableUnits.push_back("non_dim");
327  this->mInitialConditions.push_back(3.333333333333334e-03);
328 
329  this->mVariableNames.push_back("pRb_p");
330  this->mVariableUnits.push_back("non_dim");
331  this->mInitialConditions.push_back(1.000000000000000e-04);
332 
333  this->mVariableNames.push_back("APC");
334  this->mVariableUnits.push_back("non_dim");
335  this->mInitialConditions.push_back(NAN); // will be filled in later
336 
337  this->mVariableNames.push_back("Beta_Cat1");
338  this->mVariableUnits.push_back("non_dim");
339  this->mInitialConditions.push_back(NAN); // will be filled in later
340 
341  this->mVariableNames.push_back("Beta_Cat2");
342  this->mVariableUnits.push_back("non_dim");
343  this->mInitialConditions.push_back(NAN); // will be filled in later
344 
345  this->mVariableNames.push_back("Wnt");
346  this->mVariableUnits.push_back("non_dim");
347  this->mInitialConditions.push_back(NAN); // will be filled in later
348 
349  this->mInitialised = true;
350 }
351 
353 {
354  return mWntLevel;
355 }
356 
357 // Serialization for Boost >= 1.36
void SetDefaultInitialCondition(unsigned index, double initialCondition)
double CalculateRootFunction(double time, const std::vector< double > &rY)
void SetStateVariables(const std::vector< double > &rStateVariables)
void SetMutationState(boost::shared_ptr< AbstractCellMutationState > pMutationState)
bool CalculateStoppingEvent(double time, const std::vector< double > &rY)
WntCellCycleOdeSystem(double wntLevel=0.0, boost::shared_ptr< AbstractCellMutationState > pMutationState=boost::shared_ptr< AbstractCellMutationState >(), std::vector< double > stateVariables=std::vector< double >())
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
#define CHASTE_CLASS_EXPORT(T)
const boost::shared_ptr< AbstractCellMutationState > GetMutationState() const
#define NAN
Definition: IsNan.hpp:73
boost::shared_ptr< AbstractCellMutationState > mpMutationState
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)