/*
 * Mathematical model predicts a critical role for osteoclast autocrine
 * regulation in the control of bone remodeling
 * 
 * Model Status
 * 
 * This CellML model runs in both OpenCell and COR to recreate
 * the published results (the oscillating model in figure 3). The
 * parameter values have been set to those in the legend of figure
 * 3. Note there appears to be a typographical error in equation
 * A2 in the published paper: 1-(g11/gamma) has been changed to
 * (1-g11)gamma to be consistent with equations A1 and A3.
 * 
 * Model Structure
 * 
 * ABSTRACT: Bone remodeling occurs asynchronously at multiple
 * sites in the adult skeleton and involves resorption by osteoclasts,
 * followed by formation of new bone by osteoblasts. Disruptions
 * in bone remodeling contribute to the pathogenesis of disorders
 * such as osteoporosis, osteoarthritis, and Paget's disease. Interactions
 * among cells of osteoblast and osteoclast lineages are critical
 * in the regulation of bone remodeling. We constructed a mathematical
 * model of autocrine and paracrine interactions among osteoblasts
 * and osteoclasts that allowed us to calculate cell population
 * dynamics and changes in bone mass at a discrete site of bone
 * remodeling. The model predicted different modes of dynamic behavior:
 * a single remodeling cycle in response to an external stimulus,
 * a series of internally regulated cycles of bone remodeling,
 * or unstable behavior similar to pathological bone remodeling
 * in Paget's disease. Parametric analysis demonstrated that the
 * mode of dynamic behavior in the system depends strongly on the
 * regulation of osteoclasts by autocrine factors, such as transforming
 * growth factor beta. Moreover, simulations demonstrated that
 * nonlinear dynamics of the system may explain the differing effects
 * of immunosuppressants on bone remodeling in vitro and in vivo.
 * In conclusion, the mathematical model revealed that interactions
 * among osteoblasts and osteoclasts result in complex, nonlinear
 * system behavior, which cannot be deduced from studies of each
 * cell type alone. The model will be useful in future studies
 * assessing the impact of cytokines, growth factors, and potential
 * therapies on the overall process of remodeling in normal bone
 * and in pathological conditions such as osteoporosis and Paget's
 * disease.
 * 
 * model diagram
 * 
 * [[Image file: komarova_2003.png]]
 * 
 * Schematic representation of interactions between osteoclasts
 * and osteoblasts included in the model. Thick arrows represent
 * the processes of formation and removal of osteoclasts and osteoblasts.
 * Fine arrows represent the effects of autocrine and paracrine
 * regulators of bone remodeling on the rates of osteoclast and
 * osteoblast formation. TGF-beta, transforming growth factor beta,
 * released and activated by resorbing osteoclasts, directly stimulates
 * formation of osteoclasts and osteoblasts. IGF, insulin-like
 * growth factors, secreted by osteoblasts and released by resorbing
 * osteoclasts, activate osteoblast formation. RANKL, expressed
 * on and released by osteoblasts, activates osteoclastogenesis.
 * OPG, osteoprotegerin, released by osteoblasts, inhibits the
 * actions of RANKL.
 * 
 * The original paper reference is cited below:
 * 
 * Mathematical model predicts a critical role for osteoclast autocrine
 * regulation in the control of bone remodeling, Komarova SV, Smith
 * RJ, Dixon SJ, Sims SM, Wahl LM, 2003, Bone, 33, 206-215. PubMed
 * ID: 14499354
 */

import nsrunit;
// Warning: unit conversion turned off due to unit errors in 3 equation(s)
unit conversion off;
unit day=86400 second^1;
unit cell = fundamental;
//Warning:  unit percent_ renamed from percent, as the latter is predefined in JSim with different fundamental units.
unit percent_ = fundamental;
unit flux=1.1574074E-5 second^(-1)*cell^1;
unit first_order_rate_constant=1.1574074E-5 second^(-1);
unit second_order_rate_constant=1.1574074E-5 second^(-1)*cell^(-1);
unit percent_per_cell_per_day=1.1574074E-5 second^(-1)*cell^(-1)*percent_^1;

math main {
	realDomain time day;
	time.min=0;
	extern time.max;
	extern time.delta;
	real x1(time) cell;
	when(time=time.min) x1=10.06066;
	real alpha1 flux;
	alpha1=3;
	real beta1 first_order_rate_constant;
	beta1=0.2;
	real g11 dimensionless;
	g11=1.1;
	real g21 dimensionless;
	g21=-0.5;
	real x2(time) cell;
	when(time=time.min) x2=212.132;
	real alpha2 first_order_rate_constant;
	alpha2=4;
	real beta2 first_order_rate_constant;
	beta2=0.02;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real g12 dimensionless;
//	g12=1;
//	Var below replaced by constant in model eqns to satisfy unit correction
//	real g22 dimensionless;
//	g22=0;
	real z(time) percent_;
	when(time=time.min) z=100.0;
	real k1 percent_per_cell_per_day;
	k1=0.093;
	real k2 percent_per_cell_per_day;
	k2=0.0008;
	real y1(time) cell;
	real y2(time) cell;
	real x1_bar cell;
	real x2_bar cell;
	real gamma dimensionless;

	// <component name="environment">

	// <component name="x1">
	x1:time=(alpha1*x1^g11*x2^g21-beta1*x1);

	// <component name="x2">
	x2:time=(alpha2*x1^1*x2^0-beta2*x2);

	// <component name="z">
	z:time=(k2*y2-k1*y1);

	// <component name="y1">
	y1=(if (x1>x1_bar) x1-x1_bar else (0 cell));

	// <component name="y2">
	y2=(if (x2>x2_bar) x2-x2_bar else (0 cell));

	// <component name="x1_bar">
	x1_bar=((beta1/alpha1)^((1-0)/gamma)*(beta2/alpha2)^(g21/gamma));

	// <component name="x2_bar">
	x2_bar=((beta1/alpha1)^(1/gamma)*(beta2/alpha2)^((1-g11)/gamma));

	// <component name="model_parameters">
	gamma=(1*g21-(1-g11)*(1-0));
}