/* * 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; // }