/*
 * Modelling the Complex Patterns of Viral Load Decay under Antiretroviral
 * Therapy
 * 
 * Model Status
 * 
 * This model is valid CellML but can not be integrated as it is
 * overconstrained. Some variables are not initialised but more
 * importantly the model contains time delays, making it unsuitable
 * for description using CellML.
 * 
 * Model Structure
 * 
 * When treated with antiretroviral therapy, the plasma viral loads
 * in HIV patients decline in three distinct phases:
 * 
 * When the treatment is initially administered there is a 6 to
 * 72 hour delay before the viral loads begin to decline;
 * 
 * following this shoulder period viral load decline rapidly for
 * about one week;
 * 
 * finally, viral loads continue to decrease, but at a slower rate.
 * 
 * In the Dixit and Perelson 2004 publication described here, the
 * authors present a mathematical model that combines pharmacokinetics
 * and viral dynamics and includes intracellular delay. The model
 * of viral dynamics is based on that published by Perelson et
 * al. (see below). To accurately represent the pharmacokinetics
 * they assume a two compartment model consisting of blood and
 * cells. The models of viral dynamics and pharmacokinetics are
 * connected via a parameter which represents drug specific efficacy.
 * Because reverse transcriptase inhibitors (RTIs) and protease
 * inhibitors (PIs) have very different effects on the duration
 * of the intracellular delay, they have been modelled separately
 * (see and ).
 * 
 * The complete original paper reference is cited below:
 * 
 * Complex patterns of viral load decay under antiretroviral therapy:
 * influence of pharmacokinetics and intracelllar delay, Narendra
 * M. Dixit and Alan S. Perelson, 2004, Journal of Theoretical
 * Biology , 226, 95-109. PubMed ID: 14637059
 * 
 * cell diagram
 * 
 * [[Image file: dixit_2004a.png]]
 * 
 * Schematic summary of the dynamics of HIV-1 infection in vivo.
 * This model is based on that published by Perelson et al. in
 * 1996.
 * 
 * reaction diagram1
 * 
 * [[Image file: dixit_2004b.png]]
 * 
 * Pharmacokinetics - the two compartment model for protease inhibitors.
 * 
 * reaction diagram2
 * 
 * [[Image file: dixit_2004c.png]]
 * 
 * Pharmacokinetics - the two compartment model for reverse transcriptase
 * inhibitors. The main difference between the models is that RTIs
 * must be phosphorylated within the cell in order to be in their
 * active form.
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 1 equation(s)
unit conversion off;
unit ml=1E-6 meter^3;
unit mg=1E-6 kilogram^1;
unit per_ml=1E6 meter^(-3);
unit mg_per_ml=1 kilogram^1*meter^(-3);
unit day=86400 second^1;
unit flux=1.1574074E-11 meter^3*second^(-1);
unit first_order_rate_constant=1.1574074E-5 second^(-1);
unit second_order_rate_constant=11.57407407 meter^(-3)*second^(-1);

math main {
	realDomain time day;
	time.min=0;
	extern time.max;
	extern time.delta;
	real T(time) per_ml;
	when(time=time.min) T=1e6;
	real lamda second_order_rate_constant;
	lamda=1e4;
	real d first_order_rate_constant;
	d=0.01;
	real k flux;
	k=2.4e-8;
	real VI(time) per_ml;
	when(time=time.min) VI=1;
	real T_(time) per_ml;
	when(time=time.min) T_=1;
	real tau first_order_rate_constant;
	tau=1.5;
	real m first_order_rate_constant;
	m=0.01;
	real delta first_order_rate_constant;
	delta=0.01;
	real N dimensionless;
	N=2500;
	real c first_order_rate_constant;
	c=23;
	real epsilon_PI(time) dimensionless;
	real VNI(time) per_ml;
	when(time=time.min) VNI=2;
	real IC50 mg_per_ml;
	IC50=9e-7;
	real Cc(time) mg_per_ml;
	when(time=time.min) Cc=0;
	real Cb(time) mg_per_ml;
	real Vd ml;
	Vd=28000;
	real F dimensionless;
	F=1;
	real D mg;
	D=600;
	real ka first_order_rate_constant;
	ka=14.64;
	real ke first_order_rate_constant;
	ke=6.86;
	real kacell first_order_rate_constant;
	kacell=24000;
	real kecell first_order_rate_constant;
	kecell=1.1;
	real Cx(time) mg_per_ml;
	real H dimensionless;
	H=0.052;
	real fb dimensionless;
	fb=0.99;

	// <component name="environment">

	// <component name="T">
	T:time=(lamda-(d*T+k*T*VI));

	// <component name="T_">
	T_:time=(k*T*(time-tau)*VI*(time-tau)*exp((-1)*m*tau)-delta*T_);

	// <component name="VI">
	VI:time=(N*delta*T_*(1-epsilon_PI)-c*VI);

	// <component name="VNI">
	VNI:time=(N*delta*T_*epsilon_PI-c*VNI);

	// <component name="epsilon_PI">
	epsilon_PI=(Cc/(IC50+Cc));

	// <component name="Cb">
	Cb=(F*D/Vd*(ka/(ke-ka))*(exp((-1)*ka*time)-exp((-1)*ke*time)));

	// <component name="Cc">
	Cc:time=(kacell*Cx-kecell*Cc);

	// <component name="Cx">
	Cx=(if (((1-fb)*H*Cb-Cc)>(0 mg_per_ml)) (1-fb)*H*Cb-Cc else (0 mg_per_ml));

	// <component name="kinetic_parameters">
}