/*
 * A model for the dynamics of human weight cycling
 * 
 * Model Status
 * 
 * Runs in COR and OpenCell to reproduce published output (figure
 * 3). No scale is mentioned in the paper, so time is scaled arbitrarily
 * in seconds, and variables that might be represented in arbitrary
 * units of mass or force (weight) have been left dimensionless
 * as their normalisation has been mentioned in the paper. The
 * units are consistent throughout.
 * 
 * Model Structure
 * 
 * ABSTRACT: The resolution to lose weight by cognitive restraint
 * of nutritional intake often leads to repeated bouts of weight
 * loss and regain, a phenomenon known as weight cycling or "yo-yo
 * dieting". A simple mathematical model for weight cycling is
 * presented. The model is based on a feedback of psychological
 * nature by which a subject decides to reduce dietary intake once
 * a threshold weight is exceeded. The analysis of the model indicates
 * that sustained oscillations in body weight occur in a parameter
 * range bounded by critical values. Only outside this range can
 * body weight reach a stable steay state. The model provides a
 * theoretical framework that captures key facets of weight cycling
 * and suggests ways to control the phenomenon. The view that weight
 * cycling represents self-sustained oscillations has indeed specific
 * implications. In dymamical terms, to bring weight cycling to
 * an end, parameter values should change in such a way as to induced
 * the transition of body weight from sustained oscillations around
 * an unstable steady state to a stable steady state. Maintaining
 * weight under a critical value should prevent weight cycling
 * and allow body weight to stabilise below the oscillatory range.
 * 
 * A model for the dynamics of human weight cycling, Albert Goldbeter,
 * 2006, Journal of Biosciences, 31, 129-136. PubMed ID: 16595882
 * 
 * Model flowchart
 * 
 * [[Image file: goldbeter_2006.png]]
 * 
 * PQR model for weight cycling.
 */

import nsrunit;
unit conversion on;
unit s=1 second^1;
unit per_s=1 second^(-1);

math main {
	realDomain time s;
	time.min=0;
	extern time.max;
	extern time.delta;
	real P(time) dimensionless;
	when(time=time.min) P=0.0;
	real Q(time) dimensionless;
	when(time=time.min) Q=0.0;
	real a per_s;
	a=0.1;
	real b per_s;
	b=0.1;
	real K dimensionless;
	K=0.2;
	real P_star dimensionless;
	real V_3 per_s;
	V_3=6;
	real V_4 per_s;
	V_4=2.5;
	real K_3 dimensionless;
	K_3=0.01;
	real K_4 dimensionless;
	K_4=0.01;
	real R(time) dimensionless;
	when(time=time.min) R=0.0;
	real V_1 per_s;
	V_1=1;
	real V_2 per_s;
	V_2=1.5;
	real K_1 dimensionless;
	K_1=0.01;
	real K_2 dimensionless;
	K_2=0.01;
	real R_star dimensionless;

	// <component name="environment">

	// <component name="P">
	P:time=(a*Q-b*(P/(K+P)));

	// <component name="P_star">
	P_star=(V_4/V_3*((1+2*K_3)/(1+2*K_4)));

	// <component name="Q">
	Q:time=(V_1*((1-Q)/(K_1+(1-Q)))-V_2*R*(Q/(K_2+Q)));

	// <component name="R">
	R:time=(P*V_3*((1-R)/(K_3+(1-R)))-V_4*(R/(K_4+R)));

	// <component name="R_star">
	R_star=(V_1/V_2*((1+2*K_2)/(1+2*K_1)));
}