/*
 * The Follicular Automaton Model for Hair Cycles
 * 
 * Model Status
 * 
 * This CellMl model runs in OpenCell and COR to recreate the same
 * steady state output of the original published model (as seen
 * in in figure 5(d)). The units have been checked and they are
 * consistent. The CellMl model represents the deterministic version
 * of the follicular automaton model. Currently CellML cannot be
 * used to describe stochastic variables.
 * 
 * Model Structure
 * 
 * ABSTRACT: Human scalp hair consists of a set of about 10 to
 * the 5 follicles which progress independently through developmental
 * cycles. Each hair follicle successively goes through the anagen
 * (A), catagen (C), telogen (T) and latency (L) phases that correspond,
 * respectively, to growth, arrest and hair shedding before a new
 * anagen phase is initiated. Long-term experimental observations
 * in a group of ten male, alopecic and non-alopecic volunteers
 * allowed determination of the characteristics of hair follicle
 * cycles. On the basis of these observations, we previously proposed
 * a follicular automaton model to simulate the dynamics of human
 * hair cycles and the development of different patterns of alopecia
 * [Halloy et al. (2000) Proc. Natl Acad. Sci. U.S.A.97, 8328-8333].
 * The automaton model is defined by a set of rules that govern
 * the stochastic transitions of each follicle between the successive
 * states A, T, L and the subsequent return to A. These transitions
 * occur independently for each follicle, after time intervals
 * given stochastically by a distribution characterized by a mean
 * and a standard deviation. The follicular automaton model was
 * shown to account both for the dynamical transitions observed
 * in a single follicle, and for the behaviour of an ensemble of
 * independently cycling follicles. Here, we extend these results
 * and investigate additional properties of the model. We present
 * a deterministic version of the follicular automaton. We show
 * that numerical simulations of the stochastic version of the
 * automaton yield steady-state level of follicles in the different
 * phases which approach the levels predicted by the deterministic
 * equations as the number of follicles progressively increases.
 * Only the stochastic version can successfully reproduce the fluctuations
 * of the fractions of follicles in each of the three phases, observed
 * in small follicle populations. When the standard deviation is
 * reduced or when the follicles become otherwise synchronized,
 * e.g. by a periodic external signal inducing the transition of
 * anagen follicles into telogen phase, large-amplitude oscillations
 * occur in the fractions of follicles in the three phases. These
 * oscillations are not observed in humans but are reminiscent
 * of the phenomenon of moulting observed in a number of mammalian
 * species. Copyright 2002 Elsevier Science Ltd.
 * 
 * The original paper reference is cited below:
 * 
 * The Follicular Automaton Model: Effect of Stochasticity and
 * of Synchronization of Hair Cycles, J. Halloy, B.A. Bernard,
 * G. Loussouarn and A. Goldbeter, 2002, Journal of Theoretical
 * Biology , 214, 469-479. PubMed ID: 11846603
 * 
 * reaction_diagram
 * 
 * [[Image file: halloy_2002.png]]
 * 
 * The above diagram represents the transition of a model hair
 * follicle from anagen (A) to telogen (T) to latency (L) phase,
 * successively. After phase T, the follicle may either die or
 * miniaturise (transition to M; this usually occurs after a critical
 * number of cycles), or complete a cycle by entering a new A phase.
 */

import nsrunit;
unit conversion on;
unit fraction = fundamental;
unit month=2592000 second^1;
unit flux=3.8580247E-7 second^(-1)*fraction^1;
unit first_order_rate_constant=3.8580247E-7 second^(-1);

math main {
	realDomain time month;
	time.min=0;
	extern time.max;
	extern time.delta;
	real L(time) fraction;
	when(time=time.min) L=0.33;
	real mu_T month;
	mu_T=2.4;
	real mu_L month;
	mu_L=2;
	real mu_A month;
	mu_A=3.6;
	real epsilon first_order_rate_constant;
	epsilon=0;
	real M(time) fraction;
	when(time=time.min) M=0;
	real T(time) fraction;
	when(time=time.min) T=0.33;
	real L_0 dimensionless;
	real A(time) fraction;
	when(time=time.min) A=0.34;
	real A_0 dimensionless;
	real T_0 dimensionless;

	// <component name="environment">

	// <component name="L">
	L_0=(mu_L/(mu_A+mu_T+mu_L));
	L:time=(1/mu_T*T-1/mu_L*L+epsilon*M);

	// <component name="A">
	A_0=(mu_A/(mu_A+mu_T+mu_L));
	A:time=(1/mu_L*L-1/mu_A*A);

	// <component name="T">
	T_0=(mu_T/(mu_A+mu_T+mu_L));
	T:time=(1/mu_A*A-1/mu_T*T);

	// <component name="M">
	M:time=(epsilon*T);

	// <component name="reaction_constants">
}