/*
 * Influence of Delayed Viral Production on Viral Dynamics in HIV-1
 * Infected Patients
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR however we are
 * uncertain as to whether or not the CellML model replicates the
 * original model from the published paper as there are no validation
 * figures to compare it against. The CellML model is based on
 * equatiosn 10a to 10e. Parameter values for the variables k and
 * p are not stated in the paper so in the CellML model these were
 * taken from a previously published model by the same author (Perelson
 * et al. 1996). The units have been checked and they are consistent.
 * Note that while the model does run in COR, since the unit of
 * time is days the model is not ideally suited for running in
 * COR.
 * 
 * Model Structure
 * 
 * ABSTRACT: We present and analyze a model for the interaction
 * of human immunodeficiency virus type 1 (HIV-1) with target cells
 * that includes a time delay between initial infection and the
 * formation of productively infected cells. Assuming that the
 * variation among cells with respect to this 'intracellular' delay
 * can be approximated by a gamma distribution, a high flexible
 * distribution that can mimic a variety of biologically plausible
 * delays, we provide analytical solutions for the expected decline
 * in plasma virus concentration after the initiation of antiretroviral
 * therapy with one or more protease inhibitors. We then use the
 * model to investigate whether the parameters that characterize
 * viral dynamics can be identified from biological data. Using
 * non-linear least-squares regression to fit the model to simulated
 * data in which the delays conform to a gamma distribution, we
 * show that good estimates for free viral clearance rates, infected
 * cell death rates, and parameters characterizing the gamma distribution
 * can be obtained. For simulated data sets in which the delays
 * were generated using other biologically plausible distributions,
 * reasonably good estimates for viral clearance rates, infected
 * cell death rates, and mean delay times can be obtained using
 * the gamma-delay model. For simulated data sets that include
 * added simulated noise, viral clearance rate estimates are not
 * as reliable. If the mean intracellular delay is known, however,
 * we show that reasonable estimates for the viral clearance rate
 * can be obtained by taking the harmonic mean of viral clearance
 * rate estimates from a group of patients. These results demonstrate
 * that it is possible to incorporate distributed intracellular
 * delays into existing models for HIV dynamics and to use these
 * refined models to estimate the half-life of free virus from
 * data on the decline in HIV-1 RNA following treatment.
 * 
 * The original paper reference is cited below:
 * 
 * Influence of Delayed Viral Production on Viral Dynamics in HIV-1
 * Infected Patients, John E. Mittler, Bernhard Sulzer, Avidan
 * U. Neumann, and Alan S. Perelson, 1998, Mathematical Biosciences
 * , 152, 143-163. PubMed ID: 9780612
 * 
 * cell diagram
 * 
 * [[Image file: mittler_1998a.png]]
 * 
 * Schematic summary of the dynamics of HIV-1 infection in vivo
 * captured by the Perelson et al. 1996 model.
 * 
 * cell diagram
 * 
 * [[Image file: mittler_1998b.png]]
 * 
 * Schematic summary of the dynamics of viral infection in vivo
 * captured by the Herz et al. 1996 model.
 */

import nsrunit;
unit conversion on;
unit per_ml=1E6 meter^(-3);
unit day=86400 second^1;
unit first_order_rate_constant=1.1574074E-5 second^(-1);
unit ml_per_day=1.1574074E-11 meter^3*second^(-1);

math main {
	realDomain time day;
	time.min=0;
	extern time.max;
	extern time.delta;
	real T per_ml;
	real k ml_per_day;
	k=2.4e-5;
	real p first_order_rate_constant;
	p=774;
	real c first_order_rate_constant;
	c=3;
	real delta first_order_rate_constant;
	delta=0.5;
	real I(time) per_ml;
	when(time=time.min) I=0.1;
	real I_0 per_ml;
	real k_ ml_per_day;
	real E4(time) per_ml;
	when(time=time.min) E4=0;
	real VI_0 per_ml;
	VI_0=200000;
	real VI(time) per_ml;
	when(time=time.min) VI=200000;
	real h(time) dimensionless;
	real VNI(time) per_ml;
	when(time=time.min) VNI=0;
	real V(time) per_ml;
	real E1(time) per_ml;
	when(time=time.min) E1=0;
	real b_ day;
	real E2(time) per_ml;
	when(time=time.min) E2=0;
	real E3(time) per_ml;
	when(time=time.min) E3=0;
	real tau_p day;
	tau_p=0;
	real b day;
	b=0.25;
	real m first_order_rate_constant;
	m=0.01;
	real n dimensionless;
	n=4;

	// <component name="environment">

	// <component name="T">
	T=(c*delta/(k*p));

	// <component name="I">
	I:time=(k_*T*E4-delta*I);
	I_0=(c/p*VI_0);

	// <component name="VI">
	VI:time=((1-h)*p*I-c*VI);

	// <component name="VNI">
	VNI:time=(h*p*I-c*VNI);

	// <component name="virus_total">
	V=(VI+VNI);

	// <component name="E1">
	E1:time=((VI-E1)/b_);

	// <component name="E2">
	E2:time=((E1-E2)/b_);

	// <component name="E3">
	E3:time=((E2-E3)/b_);

	// <component name="E4">
	E4:time=((E3-E4)/b_);

	// <component name="Heavyside_function">
	h=(if (time<tau_p) 0 else 1);

	// <component name="kinetic_parameters">
	b_=(b/(1+m*b));
	k_=(k/(1+m*b)^n);
}