/*
* Virtual muscle: a computational approach to understanding the
* effects of muscle properties on motor control
*
* Model Status
*
* This CellML model has been checked in both COR and OpenCell.
* Currently it can be opened in COR but due to the presence of
* 'circular arguments' the model does not run. However the model
* will run in OpenCell but although it does run, it does so very,
* very slowly, and so far we have only been able to get it to
* solve for 3 iterations. The units have been checked and they
* are consistent.
*
* Model Structure
*
* Abstract: This paper describes a computational approach to modeling
* the complex mechanical properties of muscles and tendons under
* physiological conditions of recruitment and kinematics. It is
* embodied as a software package for use with Matlab™ and Simulink
* that allows the creation of realistic musculotendon elements
* for use in motor control simulations. The software employs graphic
* user interfaces (GUI) and dynamic data exchange (DDE) to facilitate
* building custom muscle model blocks and linking them to kinetic
* analyses of complete musculoskeletal systems. It is scalable
* in complexity and accuracy. The model is based on recently published
* data on muscle and tendon properties measured in feline slow-
* and fast-twitch muscle, and incorporates a novel approach to
* simulating recruitment and frequency modulation of different
* fiber-types in mixed muscles. This software is distributed freely
* over the Internet at http://ami.usc.edu/mddf/virtualmuscle.
*
* The original paper reference is cited below:
*
* Virtual muscle: a computational approach to understanding the
* effects of muscle properties on motor control, Ernest J. Cheng,
* Ian E. Brown and Gerald E. Loeb, 2000, Journal of Neuroscience
* Methods , 101, 117-130. PubMed ID: 10996372
*
* reaction diagram
*
* [[Image file: cheng_2000.png]]
*
* Schematic representation of the model equations and terms. These
* elements were designed to have a one-to-one correspondence with
* the physiological substrates of muscle contraction. The behavior
* of each element is governed by an equation driven by one to
* four input variables, with one to seven user-modifiable coefficients.
* The coefficients can be modified in the BuildFiberTypes function.
* Complete descriptions of these elements can be found in Brown
* and Loeb, 2000 and Brown. FPE1 represents the passive visco-elastic
* properties of stretching a muscle. FPE2 represents the passive
* resistance to compression of the thick filaments at short muscle
* lengths. FL represents the tetanic, isometric force–length relationship.
* FV represents the tetanic force–velocity (FV) relationship.
* Af represents the isometric, activation–frequency (Af) relationship.
* feff represents the time lag between changes in firing frequency
* and internal activation (i.e. rise and fall times). Leff represents
* the time lag between changes in length and the effect of length
* on the Af relationship. S represents the effects of ‘sag’ on
* the activation during a constant stimulus frequency. Y represents
* the effects of yielding (on activation) following movement during
* sub-maximal activation.
*/
import nsrunit;
unit conversion on;
// unit millisecond predefined
unit first_order_rate_constant=1E3 second^(-1);
math main {
realDomain time millisecond;
time.min=0;
extern time.max;
extern time.delta;
real F_SE newton;
real cT dimensionless;
cT=27.8;
real kT dimensionless;
kT=0.0047;
real LT_r dimensionless;
LT_r=0.964;
real LT dimensionless;
LT=0.02;
real F_max newton;
F_max=23;
real F_PE1(time) dimensionless;
real c1 dimensionless;
c1=23;
real k1 dimensionless;
k1=0.046;
real L_r1 dimensionless;
L_r1=1.17;
real eta millisecond;
eta=0.001;
real L(time) dimensionless;
when(time=time.min) L=0.15;
real L_max dimensionless;
L_max=0.13;
real V(time) first_order_rate_constant;
when(time=time.min) V=0.09314;
real F_PE2(time) dimensionless;
real c2 dimensionless;
c2=23;
real k2 dimensionless;
k2=0.046;
real L_r2 dimensionless;
L_r2=1.17;
real FL(time) dimensionless;
real beta dimensionless;
beta=1.55;
real omega dimensionless;
omega=0.75;
real rho dimensionless;
rho=2.12;
real FV(time) dimensionless;
real av0 dimensionless;
av0=-1.53;
real av1 dimensionless;
av1=0;
real av2 dimensionless;
av2=0;
real cv0 dimensionless;
cv0=-5.7;
real cv1 dimensionless;
cv1=9.18;
real bv first_order_rate_constant;
bv=0.69;
real V_max first_order_rate_constant;
V_max=-9.15;
real Af(time) dimensionless;
real af dimensionless;
af=0.56;
real nf0 dimensionless;
nf0=2.1;
real nf1 dimensionless;
nf1=3.3;
real nf dimensionless;
nf=1;
real Y(time) dimensionless;
when(time=time.min) Y=1;
real S(time) dimensionless;
when(time=time.min) S=1;
real f_eff(time) dimensionless;
when(time=time.min) f_eff=0;
real L_eff(time) dimensionless;
when(time=time.min) L_eff=0.1497;
real F0(time) dimensionless;
real F_CE(time) newton;
real F_total(time) newton;
real T_L millisecond;
T_L=0.088;
real T_s millisecond;
T_s=43;
real as1 dimensionless;
as1=1.76;
real as2 dimensionless;
as2=0.96;
real as_(time) dimensionless;
real c_Y dimensionless;
c_Y=0.35;
real V_Y first_order_rate_constant;
V_Y=0.1;
real T_Y millisecond;
T_Y=200;
real f_int(time) dimensionless;
when(time=time.min) f_int=0;
real df_eff_dt(time) first_order_rate_constant;
real T_f(time) millisecond;
real T_f1 millisecond;
T_f1=0.35;
real T_f2 millisecond;
T_f2=0.1;
real T_f3 millisecond;
T_f3=200;
real T_f4 millisecond;
T_f4=200;
real f_env dimensionless;
f_env=1;
real mass kilogram;
mass=0.005;
real V0(time) first_order_rate_constant;
real L0(time) dimensionless;
//
//
F_SE=(cT*F_max*kT*ln(exp((LT-LT_r)/kT)+1));
//
F_PE1=(c1*k1*ln(exp((L/L_max-L_r1)/k1)+1)+eta*V);
//
F_PE2=(c2*(exp(k2*(L-L_r2))-1));
//
FL=exp((-1)*abs((L^beta-1)/omega)^rho);
//
FV=(if (V<=(0 first_order_rate_constant)) (V_max-V)/(V_max+(cv0+cv1*L)*V) else (bv-(av0+av1*L+av2*L^2)*V)/(bv+V));
//
Af=(1-exp((-1)*(Y*S*f_eff/(af*nf))^nf));
//
F0=(Af*(FL+FV+F_PE2)+F_PE1);
//
F_CE=(F0*F_max);
//
F_total=(F_SE-F_CE);
//
L_eff:time=((L-L_eff)^3/(T_L*(1-Af)));
//
S:time=((as_-S)/T_s);
as_=(if (f_eff<.1) as1 else as2);
//
Y:time=((1-(c_Y*(1-exp((-1)*abs(V)/V_Y))+Y))/T_Y);
//
f_int:time=((f_env-f_int)/T_f);
f_eff:time=df_eff_dt;
df_eff_dt=((f_int-f_eff)/T_f);
T_f=(if (df_eff_dt>=(0 first_order_rate_constant)) T_f1*L^2+T_f2*f_env else (T_f3+T_f4*Af)/L);
//
V:time=(F_total/((1 metre)*mass));
//
V0=(V/L_max);
//
L0=(L/L_max);
//
L:time=V;
//
}