/*
 * Evolutionary dynamics of mutator phenotypes in cancer: implications
 * for chemotherapy
 * 
 * Model Status
 * 
 * This is the CellML description of the simple version of the
 * model detailed in the paper. The model runs in both OpenCell
 * and COR and the units are consistent, however it does not reproduce
 * the correct published results as some of the parameter values
 * are unknown.
 * 
 * Model Structure
 * 
 * ABSTRACT: Genetic instability is a central characteristic of
 * cancers. However, the selective forces responsible for the emergence
 * of genetic instability are not clear. We use mathematical models
 * to determine the conditions under which selection favors instability,
 * and when stable cells are advantageous. We take into account
 * the processes of DNA damage, repair, cell cycle arrest, mutation,
 * and death. We find that the rate of DNA damage can play a major
 * role in this context. In particular, an increase in the rate
 * of DNA damage can reverse the relative fitness of stable and
 * unstable cells. In terms of cancer progression, we find the
 * following results. If cells have intact apoptotic responses,
 * stable cells prevail if the DNA hit rate is low. A high DNA
 * hit rate can result in the selection of genetically unstable
 * cells. This has implications for the induction of tumors by
 * carcinogens. On the other hand, if cells are characterized by
 * impaired apoptosis, we observe the opposite. Genetic instability
 * is selected for if the DNA hit rate is low. A high DNA hit rate
 * can select against instability and result in the persistence
 * of stable cells. We propose that chemotherapy can be used to
 * reverse the relative fitness of stable and unstable cells, such
 * that unstable cells are the inferior competitors. This could
 * result in the competitive exclusion of progressing cancer cells.
 * 
 * The original paper reference is cited below:
 * 
 * Evolutionary dynamics of mutator phenotypes in cancer: implications
 * for chemotherapy, Natalia L. Komarova and Dominik Wodarz, 2003,Cancer
 * Research, 63, 6635-6642. PubMed ID: 14583456
 * 
 * cell diagram
 * 
 * [[Image file: komarova_2003.png]]
 * 
 * Schematic diagram of the model, showing the processes of cell
 * reproduction, DNA damage, repair, cell cycle arrest, mutation
 * and death.
 */

import nsrunit;
unit conversion on;
unit hour=3600 second^1;
unit first_order_rate_constant=2.7777778E-4 second^(-1);
unit per_time2=7.7160494E-8 second^(-2);
unit per_time3=2.1433471E-11 second^(-3);

math main {
	realDomain time hour;
	time.min=0;
	extern time.max;
	extern time.delta;
	real S(time) dimensionless;
	when(time=time.min) S=0.5;
	real rs first_order_rate_constant;
	rs=1;
	real epsilon_s first_order_rate_constant;
	epsilon_s=0.99;
	real alpha first_order_rate_constant;
	alpha=0.1;
	real u first_order_rate_constant;
	u=0.07;
	real beta first_order_rate_constant;
	beta=0.2;
	real phi(time) dimensionless;
	real M(time) dimensionless;
	when(time=time.min) M=0.5;
	real rm first_order_rate_constant;
	rm=1.3;
	real epsilon_m first_order_rate_constant;
	epsilon_m=0.1;
	real growth_rate first_order_rate_constant;
	real a first_order_rate_constant;
	a=0.5;

	// <component name="environment">

	// <component name="S">
	S:time=(rs*S*(1-u*(1 hour)+beta*epsilon_s*u/(1 per_time3))+alpha*u*rs*S/(1 per_time3)*((1 first_order_rate_constant)-epsilon_s)-phi*S*(1 first_order_rate_constant));

	// <component name="M">
	M:time=(rm*M*(1-u/(1 first_order_rate_constant)+beta*epsilon_m*u/(1 per_time3))+alpha*u*rm*M/(1 per_time3)*((1 first_order_rate_constant)-epsilon_m)-phi*M*(1 first_order_rate_constant));

	// <component name="phi">
	phi=(S*rs/(1 first_order_rate_constant)*(1-u/(1 first_order_rate_constant)*(1-beta*epsilon_s/(1 per_time2)-alpha/(1 first_order_rate_constant)*(1-epsilon_s/(1 first_order_rate_constant))))+M*rm/(1 first_order_rate_constant)*(1-u/(1 first_order_rate_constant)*(1-beta*epsilon_m/(1 per_time2)-alpha/(1 first_order_rate_constant)*(1-epsilon_m/(1 first_order_rate_constant)))));
	growth_rate=((rs+rm)*(1-a/(1 first_order_rate_constant)));

	// <component name="kinetic_parameters">
}