/*
 * A feedback oscillator model for circulatory autoregulation
 * 
 * Model Status
 * 
 * This CellML version of the model has been checked in COR and
 * OpenCell and the model runs to recreate the published results.
 * We are grateful to the model author for their assistance in
 * getting this CellML version of their model running. A PCEnv
 * session is also available for this model and reproduces part
 * of Figure 5 from the publication.
 * 
 * ValidateCellML confirms this model as valid CellML with full
 * unit consistency.
 * 
 * Model Structure
 * 
 * Circulatory autoregulation can be defined as the ability of
 * an isolated organ to maintain a constant, or almost constant,
 * blood flow rate over a range of perfusion pressures. In the
 * study described here Annraoi De Paor and Patrick Timmons have
 * developed a mathematical model which is based on physiological
 * data. Simulation results demonstrate that autoregulation can
 * be achieved by pressure-induced oscillations in the arteriolar
 * radius.
 * 
 * model diagram
 * 
 * [[Image file: depaor_1986.png]]
 * 
 * Schematic diagram of the De Paor et al model.
 * 
 * The complete original paper reference is cited below:
 * 
 * A feedback oscillator model for circulatory autoregulation,
 * Annraoi M. De Paor and Patrick Timmons, 1986, International
 * Journal of Control , 43, 679-688.
 */

import nsrunit;
unit conversion on;

math main {
	realDomain time dimensionless;
	time.min=0;
	extern time.max;
	extern time.delta;
	real f(time) dimensionless;
	real a dimensionless;
	a=1.0;
	real r(time) dimensionless;
	when(time=time.min) r=0.5;
	real p dimensionless;
	p=1.05;
	real beta dimensionless;
	beta=1.0;
	real alpha dimensionless;
	alpha=1.0;
	real ur0(time) dimensionless;
	real m(time) dimensionless;
	real r0 dimensionless;
	r0=0.5;
	real y(time) dimensionless;
	real r1 dimensionless;
	r1=1.2;
	real ur1(time) dimensionless;
	real z(time) dimensionless;
	real t1 dimensionless;
	t1=0.1;
	real x1(time) dimensionless;
	when(time=time.min) x1=0.0;
	real q(time) dimensionless;
	real d dimensionless;
	d=5.0;
	real k dimensionless;
	k=10.5;
	real uz(time) dimensionless;
	real x2(time) dimensionless;
	when(time=time.min) x2=0.0;
	real t2 dimensionless;
	t2=0.5;
	real x3(time) dimensionless;
	when(time=time.min) x3=0.0;
	real phi(time) dimensionless;
	real h(time) dimensionless;
	when(time=time.min) h=0.0;
	real t4 dimensionless;
	t4=20.0;

	// <component name="environment">

	// <component name="f">
	f=(a*p*r^4);

	// <component name="r">
	r:time=(beta*(p*r-(alpha*(r-r0)^2*ur0+m)));

	// <component name="y">
	y=((r-r1)*ur1);

	// <component name="z">
	z=((y-x1)/t1);

	// <component name="x1">
	x1:time=((y-x1)/t1);

	// <component name="q">
	q=(k*(1-exp((-1)*(d*z)))*uz);

	// <component name="x2">
	x2:time=((q-x2)/t2);

	// <component name="x3">
	x3:time=((x2-x3)/t2);

	// <component name="m">
	m=(x3*phi);

	// <component name="phi">
	phi=(if (r<.25) 0 else if (r>2) 0 else 1-1.3061225*(r-1.125)*(r-1.125));

	// <component name="ur0">
	ur0=(if (r>r0) 1 else 0);

	// <component name="ur1">
	ur1=(if (r>r1) 1 else 0);

	// <component name="uz">
	uz=(if (z>0) 1 else 0);

	// <component name="h">
	h:time=((f-h)/t4);

	// <component name="model_parameters">
}