% Implements rat version of myofilament model % Equations for an ODE-based myofilament model % This is frozen version for publication on 1/8/2007 % % Developed by J. Jeremy Rice, IBM Research and % Pieter de Tombe, University of Illinois, Chicago % % Developed by Jason Yang, University of Virgina, 2007 % Cardiac Systems Biology Lab % Robert M. Berne Cardiovascular Research Center % Department of Biomedical Engineering % University of Virginia % PO Box 800759, Health System % Charlottesville, VA 22908 % email: jhyang@virginia.edu % Minor modifications by Jeremy Rice % Copyright 2007 function dydt = crossbridgeODEfile(t,y_crossbridge,p_crossbridge) global timeStep Cafluxes XBforces %%%%%%%%%%%%%%%%%%%%%% % Parameters % %%%%%%%%%%%%%%%%%%%%%% % Sarcomere Geometry SLmax = 2.4; % (um) maximum sarcomere length SLmin = 1.4; % (um) minimum sarcomere length len_thin = 1.2; % (um) thin filament length len_thick = 1.65; % (um) thick filament length len_hbare = 0.1; % (um) length of bare portion of thick filament % Temperature Dependence Qkon = 1.5; Qkoff = 1.3; Qkn_p = 1.6; Qkp_n = 1.6; Qfapp = 6.25; Qgapp = 2.5; Qhf = 6.25; Qhb = 6.25; Qgxb = 6.25; % Ca binding to troponin kon = 50e-3; % (1/[ms uM]) koffL = 250e-3; % (1/ms) koffH = 25e-3; % (1/ms) perm50 = 0.5; % perm variable that controls n to p transition nperm = 15; % in Hill-like fashion kn_p = 500e-3; % (1/ms) kp_n = 50e-3; % (1/ms) koffmod = 1.0; % mod to change species % Thin filament regulation and crossbridge cycling fapp = 500e-3; % (1/ms) XB on rate gapp = 70e-3; % (1/ms) XB off rate gslmod = 6; % controls SL effect on gapp hf = 2000e-3; % (1/ms) rate between pre-force and force states hfmdc = 5; % hb = 400e-3; % (1/ms) rate between pre-force and force states hbmdc = 0; % gxb = 70e-3; % (1/ms) ATP consuming transition rate sigmap = 8; % distortion dependence of STP using transition gxb sigman = 1; % xbmodsp = 4/3; % mouse specific modification for XB cycling rates % Mean strain of strongly-bound states x_0 = 0.007; % (um) strain induced by head rotation xPsi = 2; % scaling factor balancing SL motion and XB cycling % Normalized active and passive force SLrest = 1.85; % (um) rest SL length for 0 passive force PCon_t = 0.002; % (norm Force) passive force due to titin PExp_t = 10; % these apply to trabeculae and single cells only SL_c = 2.25; % (um) resting length for collagen PCon_c = 0.02; % (norm Force) passive force due to collagen PExp_c = 70; % these apply to trabeculae and single cells only % Calculation of complete muscle response massf = 0.00005e6; % ([norm Force ms^2]/um) muscle mass visc = 0.003e3; % ([norm Force ms]/um) muscle viscosity KSE = 1; % (norm Force/um) series elastic element kxb = 120; % (mN/mm^2) maximal force Trop_conc = 70; % (uM) troponin concentration %%%%%%%%%%%%%%%%%%%%% % Variables % %%%%%%%%%%%%%%%%%%%%% % State Variables N_NoXB = y_crossbridge(1); % P_NoXB = y_crossbridge(2); % N = y_crossbridge(3); % % P = y_crossbridge(4); % XBprer = y_crossbridge(5); % XBpostr = y_crossbridge(6); % P = 1-N-XBprer-XBpostr; % SL = y_crossbridge(7); % xXBpostr= y_crossbridge(8); % xXBprer = y_crossbridge(9); % TRPNCaL = y_crossbridge(10); % Low-affinity Ca-bound troponin (uM) TRPNCaH = y_crossbridge(11); % High-affinity Ca-bound troponin (uM) intf = y_crossbridge(12); % integrated force % Time-Varying Parameters Cai = p_crossbridge(1); % intracellular Ca concentration (uM) Temp = p_crossbridge(2); % Temperature (K) contrflag = p_crossbridge(3); % 0 for isometric contraction 1 for active SLset = p_crossbridge(4); % initial SL or SL for isometric contractions SEon = p_crossbridge(5); % 0 w/o series elastic element, 1 w/element % koffmod = p_crossbridge(6); % variable parameter %%%%%%%%%%%%%%%%%%%%% % Equations % %%%%%%%%%%%%%%%%%%%%% % Compute single overlap fractions sovr_ze = min(len_thick/2,SL/2); % z-line end sovr_cle=max(SL/2-(SL-len_thin),len_hbare/2); % centerline of end len_sovr = sovr_ze-sovr_cle; % single overlap length SOVFThick = len_sovr*2/(len_thick-len_hbare); % thick filament overlap frac SOVFThin = len_sovr/len_thin; % thin filament overlap frac % Compute combined Ca binding to high- (w/XB) and low- (no XB) sites Tropreg = (1-SOVFThin)*TRPNCaL + SOVFThin*TRPNCaH; permtot = sqrt(1/(1+(perm50/Tropreg)^nperm)); inprmt = min(1/permtot, 100); % Adjustments for Ca activation, temperature, SL, stress and strain konT = kon*Qkon^((Temp-310)/10); koffLT = koffL*Qkoff^((Temp-310)/10)*koffmod; koffHT = koffH*Qkoff^((Temp-310)/10)*koffmod; kn_pT = kn_p*permtot*Qkn_p^((Temp-310)/10); kp_nT = kp_n*inprmt*Qkp_n^((Temp-310)/10); fappT = fapp*xbmodsp*Qfapp^((Temp-310)/10); gapslmd = 1 + (1-SOVFThick)*gslmod; gappT = gapp*gapslmd*xbmodsp*Qgapp^((Temp-310)/10); hfmd = exp(-sign(xXBprer)*hfmdc*((xXBprer/x_0)^2)); hbmd = exp(sign((xXBpostr-x_0))*hbmdc*(((xXBpostr-x_0)/x_0)^2)); hfT = hf*hfmd*xbmodsp*Qhf^((Temp-310)/10); hbT = hb*hbmd*xbmodsp*Qhb^((Temp-310)/10); gxbmd = heav(x_0-xXBpostr)*exp(sigmap*((x_0-xXBpostr)/x_0)^2)+... (1-heav(x_0-xXBpostr))*exp(sigman*(((xXBpostr-x_0)/x_0)^2)); gxbT = gxb*gxbmd*xbmodsp*Qgxb^((Temp-310)/10); % Regulation and corssbridge cycling state derivatives dTRPNCaL = konT*Cai*(1-TRPNCaL) - koffLT*TRPNCaL; dTRPNCaH = konT*Cai*(1-TRPNCaH) - koffHT*TRPNCaH; dN_NoXB = -kn_pT*N_NoXB + kp_nT*P_NoXB; dP_NoXB = -kp_nT*P_NoXB + kn_pT*N_NoXB; dN = -kn_pT*N + kp_nT*P; % dP = -kp_nT*P + kn_pT*N - fappT*P + gappT*XBprer + gxbT*XBpostr; dXBprer = fappT*P - gappT*XBprer - hfT*XBprer + hbT*XBpostr; dXBpostr = hfT*XBprer - hbT*XBpostr - gxbT*XBpostr; dP = -(dN+dXBprer+dXBpostr); % steady-state fractions in XBprer and XBpostr using King-Altman rule SSXBprer = (hb*fapp+gxb*fapp)/... (gxb*hf+fapp*hf+gxb*gapp+hb*fapp+hb*gapp+gxb*fapp); SSXBpostr = fapp*hf/(gxb*hf+fapp*hf+gxb*gapp+hb*fapp+hb*gapp+gxb*fapp); % normalization for scaling active and passive force (maximal force) Fnordv = kxb*x_0*SSXBpostr; % Calculate Forces (active, passive, preload, afterload) force = kxb*SOVFThick*(xXBpostr*XBpostr+xXBprer*XBprer); active = force/Fnordv; ppforce_t = sign(SL-SLrest)*PCon_t*(exp(PExp_t*abs(SL-SLrest))-1); ppforce_c = heav(SL-SL_c)*PCon_c*(exp(PExp_c*abs(SL-SL_c))-1); % ppforce_c = 0; ppforce = ppforce_t + ppforce_c; preload = sign(SLset-SLrest)*PCon_t*(exp(PExp_t*abs(SLset-SLrest))-1); afterload = 0; % either constant or due to series elastic element if (SEon == 1) afterload = KSE*(SLset-SL); % if ((SEon_LengthClamp==1)&&(index==(Npulses-1))), % afterload = 5000*(SLset-SL); % end end dintf = (-ppforce+preload-active+afterload); % total force dSL = contrflag*((intf+(SLset-SL)*visc)/massf)*... % change in SL heav(SL-SLmin)*heav(SLmax-SL); % Mean strain of strongly-bound states due to SL motion and XB cycling dutyprer = (hbT*fappT+gxbT*fappT)/... % duty fractions using the (fappT*hfT+gxbT*hfT+gxbT*gappT+hbT*fappT+hbT*gappT+gxbT*fappT); dutypostr = fappT*hfT/... % King-Alman Rule (fappT*hfT+gxbT*hfT+gxbT*gappT+hbT*fappT+hbT*gappT+gxbT*fappT); dxXBprer = dSL/2+xPsi/dutyprer*(-xXBprer*fappT+(xXBpostr-x_0-xXBprer)*hbT); dxXBpostr = dSL/2+xPsi/dutypostr*(x_0+xXBprer-xXBpostr)*hfT; % Ca buffering by low-affinity troponin C (LTRPNCa) FrSBXB = (XBpostr+XBprer)/(SSXBpostr + SSXBprer); dFrSBXB = (dXBpostr+dXBprer)/(SSXBpostr + SSXBprer); dsovr_ze = -dSL/2*heav(len_thick-SL); dsovr_cle = -dSL/2*heav((2*len_thin-SL)-len_hbare); dlen_sovr = dsovr_ze-dsovr_cle; dSOVFThin = dlen_sovr/len_thin; dSOVFThick= 2*dlen_sovr/(len_thick-len_hbare); TropTot = Trop_conc*((1-SOVFThin)*TRPNCaL + ... SOVFThin*(FrSBXB*TRPNCaH+(1-FrSBXB)*TRPNCaL)); dTropTot= Trop_conc*(-dSOVFThin*TRPNCaL+(1-SOVFThin)*dTRPNCaL + ... dSOVFThin*(FrSBXB*TRPNCaH+(1-FrSBXB)*TRPNCaL) + ... SOVFThin*(dFrSBXB*TRPNCaH+FrSBXB*dTRPNCaH-dFrSBXB*TRPNCaL+... (1-FrSBXB)*dTRPNCaL)); dforce = kxb*dSOVFThick*(xXBpostr*XBpostr+xXBprer*XBprer) + ... kxb*SOVFThick*(dxXBpostr*XBpostr+xXBpostr*dXBpostr + ... dxXBprer*XBprer+xXBprer*dXBprer); % Force Derivatives % dforce = kxb*dSOVFThick*(xXBpostr*XBpostr+xXBprer*XBprer)+... % kxb*SOVFThick*(dxXBpostr*XBpostr+xXBpostr*dXBpostr+... % dxXBprer*XBprer+xXBprer*dXBprer); dactive = 0.5*dforce/Fnordv; dppforce_t = sign(SL-SLrest)*PCon_t*PExp_t*dSL*exp(PExp_t*abs(SL-SLrest)); dppforce_c = heav(SL-SL_c)*PCon_c*PExp_c*dSL*exp(PExp_c*abs(SL-SL_c)); % ppforce_c = 0; dppforce = dppforce_t + dppforce_c; dsfib = dppforce+dactive; % Update state derivatives dydt = zeros(size(y_crossbridge)); dydt(1) = dN_NoXB; dydt(2) = dP_NoXB; dydt(3) = dN; dydt(4) = dP; dydt(5) = dXBprer; dydt(6) = dXBpostr; dydt(7) = dSL; dydt(8) = dxXBpostr; dydt(9) = dxXBprer; dydt(10) = dTRPNCaL; dydt(11) = dTRPNCaH; dydt(12) = dintf; dydt(12) = 0; Cafluxes(timeStep,4:6) = [dTropTot TropTot Tropreg]; XBforces(timeStep,1:7) = [ppforce preload active afterload -dintf ... (ppforce+active) dsfib]; % XBforces(timeStep,1:7) = [ppforce_t ppforce_c dppforce_t dppforce_c -dintf ... % (ppforce+active) dppforce];