/*
 * A Mathematical Model on Germinal Center Kinetics and Termination
 * 
 * Model Status
 * 
 * This model runs in COR and OpenCell with no errors. The units
 * have been checked and they are consistent. The CellML model
 * can replicate the shapes of the graphs of all of the figures
 * in the paper but the scale does not always match those in the
 * paper.
 * 
 * Model Structure
 * 
 * ABSTRACT: We devise a mathematical model to study germinal center
 * (GC) kinetics. Earlier models for GC kinetics are extended by
 * explicitly modeling 1) the cell division history of centroblasts,
 * 2) the Ag uptake by centrocytes, and 3) T cell dynamics. Allowing
 * for T cell kinetics and T-B cell interactions, we study the
 * role of GC T cells in GC kinetics, GC termination, and B cell
 * selection. We find that GC T cells play a major role in GC formation,
 * but that the maintenance of established GC reactions requires
 * very few T cells only. The results therefore suggest that the
 * termination of a GC reaction is largely caused by lack of Ag
 * on the follicular dendritic cells and is hardly influenced by
 * Th cells. Ag consumption by centrocytes is the major factor
 * determining the decay rate of the antigenic stimulus during
 * a GC reaction. Investigating the effect of the Ag dose on GC
 * kinetics, we find that both the total size of the GC and its
 * duration are hardly influenced by the initial amount of Ag.
 * In the model this is due to a buffering effect by competition
 * for limited T cell help and/or competition between proliferating
 * centroblasts.
 * 
 * The original paper reference is cited below:
 * 
 * A Mathematical Model on Germinal Center Kinetics and Termination,
 * Can Kesmir and Rob J. De Boer, 1999, The Journal of Immunology,
 * 163, 2463-2469. PubMed ID: 10452981
 * 
 * cell diagram
 * 
 * [[Image file: kesmir_1999.png]]
 * 
 * The model GC reaction.
 */

import nsrunit;
unit conversion on;
unit cells_per_GC = fundamental;
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 B0(time) cells_per_GC;
	when(time=time.min) B0=3;
	real pr dimensionless;
	pr=0.15;
	real mu first_order_rate_constant;
	mu=3;
	real rho first_order_rate_constant;
	rho=4;
	real delta_B first_order_rate_constant;
	delta_B=0.8;
	real CT_star(time) cells_per_GC;
	real B1(time) cells_per_GC;
	when(time=time.min) B1=0;
	real alpha_B(time) dimensionless;
	real B2(time) cells_per_GC;
	when(time=time.min) B2=0;
	real B3(time) cells_per_GC;
	when(time=time.min) B3=0;
	real B4(time) cells_per_GC;
	when(time=time.min) B4=0;
	real B5(time) cells_per_GC;
	when(time=time.min) B5=0;
	real B6(time) cells_per_GC;
	when(time=time.min) B6=0;
	real B7(time) cells_per_GC;
	when(time=time.min) B7=0;
	real B8(time) cells_per_GC;
	when(time=time.min) B8=0;
	real B9(time) cells_per_GC;
	when(time=time.min) B9=0;
	real B10(time) cells_per_GC;
	when(time=time.min) B10=0;
	real B_sum(time) cells_per_GC;
	real C(time) cells_per_GC;
	when(time=time.min) C=0;
	real d first_order_rate_constant;
	d=2;
	real C_star(time) cells_per_GC;
	when(time=time.min) C_star=0;
	real CA(time) cells_per_GC;
	real C_starsum(time) cells_per_GC;
	real M(time) cells_per_GC;
	when(time=time.min) M=0;
	real A(time) cells_per_GC;
	when(time=time.min) A=500;
	real z first_order_rate_constant;
	z=0.02;
	real u dimensionless;
	u=0.15;
	real log_A(time) dimensionless;
	real T(time) cells_per_GC;
	when(time=time.min) T=0;
	real p first_order_rate_constant;
	p=2;
	real sigma first_order_rate_constant;
	sigma=5;
	real delta_T first_order_rate_constant;
	delta_T=0.8;
	real alpha_T(time) dimensionless;
	real SA dimensionless;
	SA=500;
	real ST dimensionless;
	ST=50;
	real KB dimensionless;
	KB=1e4;
	real KT dimensionless;
	KT=100;
	real total(time) cells_per_GC;
	real log_total(time) dimensionless;

	// <component name="environment">

	// <component name="B0">
	B0:time=(pr*mu*CT_star-(rho*B0+delta_B*B0));

	// <component name="B1">
	B1:time=(rho*(1+alpha_B)*B0-(rho*B1+delta_B*B1));

	// <component name="B2">
	B2:time=(rho*(1+alpha_B)*B1-(rho*B2+delta_B*B2));

	// <component name="B3">
	B3:time=(rho*(1+alpha_B)*B2-(rho*B3+delta_B*B3));

	// <component name="B4">
	B4:time=(rho*(1+alpha_B)*B3-(rho*B4+delta_B*B4));

	// <component name="B5">
	B5:time=(rho*(1+alpha_B)*B4-(rho*B5+delta_B*B5));

	// <component name="B6">
	B6:time=(rho*(1+alpha_B)*B5-(rho*B6+delta_B*B6));

	// <component name="B7">
	B7:time=(rho*(1+alpha_B)*B6-(rho*B7+delta_B*B7));

	// <component name="B8">
	B8:time=(rho*(1+alpha_B)*B7-(rho*B8+delta_B*B8));

	// <component name="B9">
	B9:time=(rho*(1+alpha_B)*B8-(rho*B9+delta_B*B9));

	// <component name="B10">
	B10:time=(rho*(1+alpha_B)*B9-(rho*B10+delta_B*B10));

	// <component name="centroblasts_sum">
	B_sum=(B1+B1+B3+B4+B5+B6+B6+B7+B8+B9+B10);

	// <component name="C">
	C:time=(d*rho*B10*(1 day)-mu*C);

	// <component name="C_star">
	C_star:time=(mu*CA-mu*C_star);

	// <component name="centrocytes_sum">
	C_starsum=(C+C_star);

	// <component name="M">
	M:time=((1-pr)*mu*CT_star);

	// <component name="A">
	A:time=((-1)*z*A-u*CA*(1 first_order_rate_constant));
	log_A=log(A/(1 cells_per_GC));

	// <component name="T">
	T:time=(sigma*(1 cells_per_GC)+p*alpha_T*CT_star-delta_T*T);

	// <component name="CA">
	CA=(C*A/(SA*(1 cells_per_GC)+A));

	// <component name="CT_star">
	CT_star=(C_star*T/(ST*(1 cells_per_GC)+C_star));

	// <component name="alpha_B">
	alpha_B=(KB/(KB+B_sum/(1 cells_per_GC)));

	// <component name="alpha_T">
	alpha_T=(KT/(KT+T/(1 cells_per_GC)));

	// <component name="total">
	total=(B_sum+C_starsum);
	log_total=log(total/(1 cells_per_GC)+1E-12);

	// <component name="kinetic_parameters">
}