/*
 * Self-tolerance and Autoimmunity in a Regulatory T Cell Model
 * 
 * Model Status
 * 
 * This CellML model represents system 4 (equations 15a-15d) in
 * the original publication. The model runs in both COR and OpenCell
 * to replicate the published results (figure 2) with R0=0.8. The
 * units have been checked and they are consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: The class of immunosuppressive lymphocytes known as
 * regulatory T cells (Tregs) has been identified as a key component
 * in preventing autoimmune diseases. Although Tregs have been
 * incorporated previously in mathematical models of autoimmunity,
 * we take a novel approach which emphasizes the importance of
 * professional antigen presenting cells (pAPCs). We examine three
 * possible mechanisms of Treg action (each in isolation) through
 * ordinary differential equation (ODE) models. The immune response
 * against a particular autoantigen is suppressed both by Tregs
 * specific for that antigen and by Tregs of arbitrary specificities,
 * through their action on either maturing or already mature pAPCs
 * or on autoreactive effector T cells. In this deterministic approach,
 * we find that qualitative long-term behaviour is predicted by
 * the basic reproductive ratio R(0) for each system. When R(0)
 * is less tHAN 1, only the trivial equilibrium exists and is stable;
 * when R(0) is greater than 1, this equilibrium loses its stability
 * and a stable non-trivial equilibrium appears. We interpret the
 * absence of self-damaging populations at the trivial equilibrium
 * to imply a state of self-tolerance, and their presence at the
 * non-trivial equilibrium to imply a state of chronic autoimmunity.
 * Irrespective of mechanism, our model predicts that Tregs specific
 * for the autoantigen in question play no role in the system's
 * qualitative long-term behaviour, but have quantitative effects
 * that could potentially reduce an autoimmune response to sub-clinical
 * levels. Our results also suggest an important role for Tregs
 * of arbitrary specificities in modulating the qualitative outcome.
 * A stochastic treatment of the same model demonstrates that the
 * probability of developing a chronic autoimmune response increases
 * with the initial exposure to self antigen or autoreactive effector
 * T cells. The three different mechanisms we consider, while leading
 * to a number of similar predictions, also exhibit key differences
 * in both transient dynamics (ODE approach) and the probability
 * of chronic autoimmunity (stochastic approach).
 * 
 * model diagram
 * 
 * [[Image file: alexander_2010.png]]
 * 
 * Flow chart illustrating interactions among populations and flow
 * in/out of compartments in System 1. Populations to be modelled
 * explicitly are in black boxes; those that are considered only
 * implicitly or as intermediaries are in grey boxes. Movement
 * in or out of a compartment is indicated by a black arrow; activating
 * influences are indicated with green arrows; suppressive influences
 * are indicated with red arrows; and other interactions or effects
 * are indicated with blue arrows.
 * 
 * The original paper reference is cited below:
 * 
 * Self-tolerance and Autoimmunity in a Regulatory T Cell Model,
 * Alexander HK and Wahl LM, 2010, Bulletin of Mathematical Biology.
 * PubMed ID: 20195912
 */

import nsrunit;
unit conversion on;
unit day=86400 second^1;
unit first_order_rate_constant=1.1574074E-5 second^(-1);

math main {
	realDomain time day;
	time.min=0;
	extern time.max;
	extern time.delta;
	real A(time) dimensionless;
	when(time=time.min) A=1.0;
	real v_tilday first_order_rate_constant;
	v_tilday=0.0025;
	real f dimensionless;
	f=1e-4;
	real muA first_order_rate_constant;
	muA=0.25;
	real b4 dimensionless;
	b4=1.00;
	real sigma4 dimensionless;
	sigma4=1.2e-5;
	real G(time) dimensionless;
	when(time=time.min) G=1e8;
	real R(time) dimensionless;
	when(time=time.min) R=0.0;
	real pi4 first_order_rate_constant;
	pi4=0.016;
	real beta first_order_rate_constant;
	beta=200.0;
	real muR first_order_rate_constant;
	muR=0.25;
	real E(time) dimensionless;
	when(time=time.min) E=0.0;
	real lambdaE first_order_rate_constant;
	lambdaE=1000.0;
	real muE first_order_rate_constant;
	muE=0.25;
	real gamma first_order_rate_constant;
	gamma=2000.0;
	real muG first_order_rate_constant;
	muG=5.0;

	// <component name="environment">

	// <component name="A">
	A:time=(f*v_tilday*G/(1+b4+sigma4*R)-muA*A);

	// <component name="R">
	R:time=((pi4*E+beta)*A-muR*R);

	// <component name="E">
	E:time=(lambdaE*A-muE*E);

	// <component name="G">
	G:time=(gamma*E-(v_tilday*G+muG*G));
}