/*
 * Teaching from classic papers: Hill's model of muscle contraction
 * 
 * Model Status
 * 
 * This CellML model runs in both PCEnv and COR to reproduce the
 * published results. The units have been checked and they are
 * consistent.
 * 
 * Model Structure
 * 
 * ABSTRACT: A. V. Hill's 1938 paper "The heat of shortening and
 * the dynamic constants of muscle" is an enduring classic, presenting
 * detailed methods, meticulous experiments, and the model of muscle
 * contraction that now bears Hill's name. Pairing a simulation
 * based on Hill's model with a reading of his paper allows students
 * to follow his thought process to discover key principles of
 * muscle physiology and gain insight into how to develop quantitative
 * models of physiological processes. In this article, the experience
 * of the author using this approach in a graduate biomedical engineering
 * course is outlined, along with suggestions for adapting this
 * approach to other audiences.
 * 
 * Teaching from classic papers: Hill's model of muscle contraction,
 * J.W. Holmes, 2006, Advances in Physiology Education, 30, 67-72.
 * PubMed ID: 16709736
 * 
 * diagram
 * 
 * [[Image file: holmes_2006.png]]
 * 
 * A schematic diagram of the two component muscle model. Force
 * in the muscle = F_SE = F_CE. By specifying Lm, Fm is calculated
 * from the force velocity properties of CE and the force-length
 * relationship of SE.
 */

import nsrunit;
unit conversion on;
unit ms=.001 second^1;
unit mNpermmsq=1E3 kilogram^1*meter^(-1)*second^(-2);
unit pms=1E3 second^(-1);

math main {
	realDomain time second;
	time.min=0;
	extern time.max;
	extern time.delta;
	real a mNpermmsq;
	a=37.24;
	real b pms;
	b=0.325;
	real Po mNpermmsq;
	Po=144.9;
	real alpha mNpermmsq;
	alpha=1449.027;
	real L_se_o dimensionless;
	L_se_o=0.3;
	real L(time) dimensionless;
	real v(time) pms;
	real L_se(time) dimensionless;
	real L_ce(time) dimensionless;
	when(time=time.min) L_ce=0.7;
	real P(time) mNpermmsq;

	// <component name="environment">

	// <component name="contraction">
	v=((-1)*b*(Po-P)/(P+a));
	L_se=(L-L_ce);
	L=(if (time<=(1 second)) 1 else if ((time>(1 second)) and (time<(5 second))) .92 else .9);
	L_ce:time=v;
	P=(alpha*(L_se-L_se_o));
}