/*
 * Parathyroid hormone temporal effects on bone formation and resorption
 * 
 * Model Status
 * 
 * This CellML model runs in OpenCell to recreate the published
 * results. The units have been checked and they are consistent.
 * PTH is pulsed every 6 hours and the model needs to be run for
 * some time before the oscillations stabilise. The CellML model
 * will also run in COR but due to the time unit being hours it
 * is not really suitable for full length simulations.
 * 
 * Model Structure
 * 
 * ABSTRACT: Parathyroid hormone (PTH) paradoxically causes net
 * bone loss (resorption) when administered in a continuous fashion,
 * and net bone formation (deposition) when administered intermittently.
 * Currently no pharmacological formulations are available to promote
 * bone formation, as needed for the treatment of osteoporosis.
 * The paradoxical behavior of PTH confuses endocrinologists, thus,
 * a model bone resorption or deposition dependent on the timing
 * of PTH administration would de-mystify this behavior and provide
 * the basis for logical drug formulation. We developed a mathematical
 * model that accounts for net bone loss with continuous PTH administration
 * and net bone formation with intermittent PTH administration,
 * based on the differential effects of PTH on the osteoblastic
 * and osteoclastic populations of cells. Bone, being a major reservoir
 * of body calcium, is under the hormonal control of PTH. The overall
 * effect of PTH is to raise plasma levels of calcium, partly through
 * bone resorption. Osteoclasts resorb bone and liberate calcium,
 * but they lack receptors for PTH. The preosteoblastic precursors
 * and preosteoblasts possess receptors for PTH, upon which the
 * hormone induces differentiation from the precursor to preosteoblast
 * and from the preosteoblast to the osteoblast. The osteoblasts
 * generate IL-6; IL-6 stimulates preosteoclasts to differentiate
 * into osteoclasts. We developed a mathematical model for the
 * differentiation of osteoblastic and osteoclastic populations
 * in bone, using a delay time of 1 hour for differentiation of
 * preosteoblastic precursors into preosteoblasts and 2 hours for
 * the differentiation of preosteoblasts into osteoblasts. The
 * ratio of the number of osteoblasts to osteoclasts indicates
 * the net effect of PTH on bone resorption and deposition; the
 * timing of events producing the maximum ratio would induce net
 * bone deposition. When PTH is pulsed with a frequency of every
 * hour, the preosteoblastic population rises and decreases in
 * nearly a symmetric pattern, with 3.9 peaks every 24 hours, and
 * 4.0 peaks every 24 hours when PTH is administered every 6 hours.
 * Thus, the preosteoblast and osteoblast frequency depends more
 * on the nearly constant value of the PTH, rather than on the
 * frequency of the PTH pulsations. Increasing the time delay gradually
 * increases the mean value for the number of osteoblasts. The
 * osteoblastic population oscillates for all intermittent administrations
 * of PTH and even when the PTH infusion is constant. The maximum
 * ratio of osteoblasts to osteoclasts occurs when PTH is administered
 * in pulses of every 6 hours. The delay features in the model
 * bear most of the responsibility for the occurrence of these
 * oscillations, because without the delay and in the presence
 * of constant PTH infusions, no oscillations occur. However, with
 * a delay, under constant PTH infusions, the model generates oscillations.
 * The osteoblast oscillations express limit cycle behavior. Phase
 * plane analysis show simple and complex attractors. Subsequent
 * to a disturbance in the number of osteoblasts, the osteoblasts
 * quickly regain their oscillatory behavior and cycle back to
 * the original attractor, typical of limit cycle behavior. Further,
 * because the model was constructed with dissipative and nonlinear
 * features, one would expect ensuing oscillations to show limit
 * cycle behavior. The results from our model, increased bone deposition
 * with intermittent PTH administration and increased bone resorption
 * with constant PTH administration, conforms with experimental
 * observations and with an accepted explanation for osteoporosis.
 * 
 * model diagram
 * 
 * [[Image file: kroll_2000.png]]
 * 
 * Schematic diagram of the effect of PTH on the development of
 * osteoblasts. PTH binds to receptors on the preosteoblast precursors
 * and stimulates their transition to preosteoblasts. However,
 * PTH binding to these preosteoblasts inhibits their differentiation
 * into osteoblasts, and IL-6 (which is secreted by the osteoblasts)
 * is believed to enhance this inhibitory effect. The osteoblasts
 * then differentiate into osteocytes at a rate which is dependent
 * on their number. The IL-6 produced by the osteoblasts also stimulates
 * the differentiation of preosteoclasts to osteoclasts. Osteoclasts
 * become senescent at a rate dependent on their number.
 * 
 * The original paper reference is cited below:
 * 
 * Parathyroid hormone temporal effects on bone formation and resorption,
 * Martin H. Kroll, 2000, Bulletin of Mathematical Biology, 62,
 * 163-187. PubMed ID: 10824426
 */

import nsrunit;
unit conversion on;
unit hour=3600 second^1;
unit first_order_rate_constant=2.7777778E-4 second^(-1);

math main {
	realDomain time hour;
	time.min=0;
	extern time.max;
	extern time.delta;
	real Y(time) dimensionless;
	when(time=time.min) Y=10.0;
	real C1 dimensionless;
	C1=1.0;
	real k1 first_order_rate_constant;
	k1=1.0;
	real P(time) dimensionless;
	when(time=time.min) P=10.0;
	real K dimensionless;
	K=5.0;
	real C dimensionless;
	C=50.0;
	real C2 dimensionless;
	C2=1.0;
	real k2 first_order_rate_constant;
	k2=1.3;
	real ky first_order_rate_constant;
	ky=0.01;
	real X(time) dimensionless;
	when(time=time.min) X=500.0;
	real k3 first_order_rate_constant;
	k3=0.05;
	real k4 first_order_rate_constant;
	k4=0.9;
	real Z(time) dimensionless;
	when(time=time.min) Z=200.0;
	real C3 dimensionless;
	C3=1.0;
	real k5 first_order_rate_constant;
	k5=5.0;
	real K2 dimensionless;
	K2=2.0;
	real IL6(time) dimensionless;
	when(time=time.min) IL6=1.9;
	real k6 first_order_rate_constant;
	k6=0.02;

	// <component name="environment">

	// <component name="Y">
	Y:time=(k1*C1*(P/(K+P))*C-(k2*C2*(1-P/(K+P))*Y+ky*Y));

	// <component name="X">
	X:time=(k2*C2*(1-P/(K+P))*Y-k3*X);

	// <component name="P">
	P:time=(if ((time>=(0 hour)) and (time<(6 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(6 hour)) and (time<(12 hour))) (-1)*(k4*P) else if ((time>=(12 hour)) and (time<(18 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(18 hour)) and (time<(24 hour))) (-1)*(k4*P) else if ((time>=(24 hour)) and (time<(30 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(30 hour)) and (time<(36 hour))) (-1)*(k4*P) else if ((time>=(36 hour)) and (time<(42 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(42 hour)) and (time<(48 hour))) (-1)*(k4*P) else if ((time>=(48 hour)) and (time<(54 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(54 hour)) and (time<(60 hour))) (-1)*(k4*P) else if ((time>=(60 hour)) and (time<(66 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(66 hour)) and (time<(72 hour))) (-1)*(k4*P) else if ((time>=(72 hour)) and (time<(78 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(78 hour)) and (time<(84 hour))) (-1)*(k4*P) else if ((time>=(84 hour)) and (time<(90 hour))) (10 first_order_rate_constant)-k4*P else if ((time>=(90 hour)) and (time<(96 hour))) (-1)*(k4*P) else 0);

	// <component name="Z">
	Z:time=(k5*C3*(IL6/(K2+IL6))-k6*Z);

	// <component name="IL6">
	IL6:time=((.1 first_order_rate_constant)*X-(10 first_order_rate_constant)*IL6);

	// <component name="model_parameters">
}