Code used to calculate time-evolution of fibril length distribution

Note this code is more general than that used to predict the cycle test data as it includes the possibility of having a non-zero starting time and/or having a gap between the end of the first cycle and the start of the second

Critical times (tA-tEE), maximum strain (emax) and rate of strain (ν)

tA=0; (*time of start of first cycle - not given a name in paper*)
tC=1/2;  (*time of peak of first cycle - called tB in paper*)
tE=1; (*time of end of first cycle - called tD in paper*)
tAA=1; (*time of start of second cycle - assumed to be same time as end of first cycle in paper*)
tCC=3/2; (*time of peak of second cycle - called tB' in paper*)
tEE=2; (*time of end of second cycle - called tD' in paper*) 
emax=0.06;(*peak fascicle strain*)
ν=emax/tC;(*strain rate*)

(*imposed fascicle strain*)
ee1[t_]:=ν t
ee2[t_]:=ν tC-ν(t-tC)
ee3[t_]:=ν(t-tAA)
ee4[t_]:=ν(tCC-tAA)-ν(t-tCC)
e[t_]:=If[t<tC,ee1[t],If[t<tE,ee2[t],If[t<tAA,0,If[t<tCC,ee3[t],If[t<tEE,ee4[t],0]]]]]

Plot[e[t],{t,0,2}]

Critical time (called tA in paper) at which fibril becomes taut

tB[ec1_]:=ec1/ν

Instantaneous (E0) and long-time (E∞) Young's modulus and relaxation time (τr)

E0=3;
E∞=1;
τr=0.3; (*assumed relaxation function parameters - for simplicity, we assume a one-term Prony series for the purposes of this code*)

Relaxation function and its derivative

EE[t_]:=E∞+(E0-E∞) E^(-t/τr)

DE[t_]:=Evaluate[D[EE[t],t]]

Plot[EE[t],{t,0,2}]

Instantaneous (J0) and long-time (J∞) creep consants and creep time (τc)

J0=1/E0;
J∞=1/E∞;
τc=τr E0/E∞;

Creep function and its derivative

J[t_]:=J∞+(J0-J∞)E^(-t/τc)

DJ[t_]:=Evaluate[D[J[t],t]]

Plot[J[t],{t,0,2}]

Expressions for σ1 (these and all following expressions were determined by analytically evaluating relevant integrals in Mathematica)

σ11[t_,ec1_]:=(((E0-E∞) (ec1-t ν+(E^((-t+tB[ec1])/(τr))) (-ec1+ν (tB[ec1]-τr))+ν τr))/(1+ec1))

σ12[t_,ec1_]:=-(((E^(-(t/(τr)))) (E0-E∞) ((E^(tB[ec1]/τr)) (ec1+ν (-tB[ec1]+τr))-(E^(tC/τr)) (ec1+ν (-tC+τr))))/(1+ec1))

σ13[t_,ec1_]:=(((E0-E∞) (ec1+t ν-2 tC ν-ν τr+(E^((-t+tC)/(τr))) (-ec1+ν (tC+τr))))/(1+ec1))

Find tD

σ1[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,EE[0]/(1+ec1)(ee1[t]-ec1)+σ11[t,ec1],EE[0]/(1+ec1)(ee2[t]-ec1)+σ12[t,ec1]+σ13[t,ec1]]]

tD[ec1_]:=tD/.FindRoot[EE[0]/(1+ec1)(ee2[tD]-ec1)+σ12[tD,ec1]+σ13[tD,ec1]==0,{tD,(tC+tE)/2}] (*called tC in paper*)

Expressions for e1

e11[t_,ec1_]:=-(1/((1+ec1) τc))(J0-J∞) ((((E^((-t+tB[ec1])/(τc))) τc (ec1+ν (-tB[ec1]+τc)) (E∞ τc-E0 τr))/(τc-τr))+(E^((-t+tC)/(τc))) τc (-ec1 E∞+E∞ tC ν-E∞ ν τc+E0 ν τr-E∞ ν τr+(((E^((tB[ec1]-tC)/(τr))) (E0-E∞) τr (ec1+ν (-tB[ec1]+τr)))/(τc-τr))))

e12[t_,ec1_]:=(1/((1+ec1) (τc-τr)))(E^(-(t/(τc)))) (J0-J∞) (-(E^(tC ((1/(τc))-(1/(τr))))) ((E0-E∞) τr (2 (E^(tC/τr)) ν τr-(E^(tB[ec1]/τr)) (ec1+ν (-tB[ec1]+τr)))+(E^(tC/τr)) (-τc+τr) (-ec1 E∞+ν (-E0 τr+E∞ (tC+τc+τr))))+(E^(tD[ec1] ((1/(τc))-(1/(τr))))) ((E0-E∞) τr (2 (E^(tC/τr)) ν τr-(E^(tB[ec1]/τr)) (ec1+ν (-tB[ec1]+τr)))+(E^(tD[ec1]/τr)) (-τc+τr) (-ec1 E∞+ν (-E0 τr+E∞ (2 tC-tD[ec1]+τc+τr)))))

Find tBB (called tA' in paper)

e1[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,1/(1+ec1)(ee1[t]-ec1),If[t<tD[ec1],1/(1+ec1)(ee2[t]-ec1),e11[t,ec1]+e12[t,ec1]]]]

tBB[ec1_]:=tBB/.FindRoot[(1+ec1)(e11[tBB,ec1]+e12[tBB,ec1])+ec1-ee3[tBB]==0,{tBB,tE}]

Find tDD (called tC' in paper)

σ2[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,EE[0]/(1+ec1)(e[t]-ec1)+σ11[t,ec1],If[t<tD[ec1],EE[0]/(1+ec1)(e[t]-ec1)+σ12[t,ec1]+σ13[t,ec1],0]]]

e2[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,1/(1+ec1)(ee1[t]-ec1),If[t<tD[ec1],1/(1+ec1)(ee2[t]-ec1),If[t<tBB[ec1],e11[t,ec1]+e12[t,ec1],If[t<tCC,1/(1+ec1)(ee3[t]-ec1),1/(1+ec1)(ee4[t]-ec1)]]]]]

Expressions for σ3

σ31[t_,ec1_]:=σ12[t,ec1]

σ32[t_,ec1_]:=(((E^(-(t/(τr)))) (E0-E∞) ((E^(tD[ec1]/τr)) (ec1+ν (-2 tC+tD[ec1]-τr))+(E^(tC/τr)) (-ec1+ν (tC+τr))))/(1+ec1))

σ33[t_,ec1_]:=(1/((1+ec1) (((τc-τr))^2)))(E^(-((t τc+tC τc+tD[ec1] τc+tBB[ec1] τr+tD[ec1] τr)/(τc τr)))) ((E^(((tD[ec1])/(τc))+((tBB[ec1])/(τr))))-(E^(((tBB[ec1])/(τc))+((tD[ec1])/(τr))))) (E0-E∞) (J0-J∞) τc (-2 (E^(((tD[ec1])/(τc))+((2 tC)/(τr)))) (E0-E∞) ν ((τr)^2)-2 (E^(tC ((1/(τc))+(1/(τr)))+((tD[ec1])/(τr)))) ν τc (E∞ τc-E0 τr)+(E^(((tB[ec1])/(τc))+(tC/(τr))+((tD[ec1])/(τr)))) (ec1+ν (-tB[ec1]+τc)) (E∞ τc-E0 τr)+(E^(((tD[ec1])/(τc))+((tB[ec1])/(τr))+(tC/(τr)))) (E0-E∞) τr (ec1+ν (-tB[ec1]+τr))-(E^((tC τc+tD[ec1] (τc+τr))/(τc τr))) (-τc+τr) (-ec1 E∞+ν (-E0 τr+E∞ (2 tC-tD[ec1]+τc+τr))))

σ34[t_,ec1_]:=(((E0-E∞) (τr (ec1+ν (-t+tAA+τr))-(E^((-t+tBB[ec1])/(τr))) τr (ec1+ν (tAA-tBB[ec1]+τr))))/((1+ec1) τr))

σ35[t_,ec1_]:=-(((E^(-(t/(τr)))) (E0-E∞) ((E^(tBB[ec1]/τr)) (ec1+ν (tAA-tBB[ec1]+τr))-(E^(tCC/τr)) (ec1+ν (tAA-tCC+τr))))/(1+ec1))

σ36[t_,ec1_]:=(((E0-E∞) ((ec1+ν (t+tAA-2 tCC-τr)) τr-(E^((-t+tCC)/(τr))) (ec1+ν (tAA-tCC-τr)) τr))/((1+ec1) τr))

σ3[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,EE[0]/(1+ec1)(e[t]-ec1)+σ11[t,ec1],If[t<tD[ec1],EE[0]/(1+ec1)(e[t]-ec1)+σ12[t,ec1]+σ13[t,ec1],If[t<tBB[ec1],0,If[t<tCC,EE[0]/(1+ec1)(ee3[t]-ec1)+σ31[t,ec1]+σ32[t,ec1]+σ33[t,ec1]+σ34[t,ec1],EE[0]/(1+ec1)(ee4[t]-ec1)+σ31[t,ec1]+σ32[t,ec1]+σ33[t,ec1]+σ35[t,ec1]+σ36[t,ec1]]]]]]

tDD[ec1_]:=tDD/.FindRoot[EE[0]/(1+ec1)(ee4[tDD]-ec1)+σ31[tDD,ec1]+σ32[tDD,ec1]+σ33[tDD,ec1]+σ35[tDD,ec1]+σ36[tDD,ec1]==0,{tDD,tCC}]

Fibril stress

σf[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,EE[0]/(1+ec1)(e[t]-ec1)+σ11[t,ec1],If[t<tD[ec1],EE[0]/(1+ec1)(e[t]-ec1)+σ12[t,ec1]+σ13[t,ec1],If[t<tBB[ec1],0,If[t<tCC,EE[0]/(1+ec1)(ee3[t]-ec1)+σ31[t,ec1]+σ32[t,ec1]+σ33[t,ec1]+σ34[t,ec1],If[t<tDD[ec1],EE[0]/(1+ec1)(ee4[t]-ec1)+σ31[t,ec1]+σ32[t,ec1]+σ33[t,ec1]+σ35[t,ec1]+σ36[t,ec1],0]]]]]]

Expressions for ef

ef1[t_,ec1_]:=(1/((1+ec1) (((τc-τr))^3)))(E^(-((tC τc+tCC τc+tD[ec1] τc+t τr+tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr)))) (J0-J∞) (-2 (E^((2 tBB[ec1] τc+2 tC τc+tCC τr+2 tD[ec1] τr)/(τc τr))) (((E0-E∞))^2) (J0-J∞) ν τc ((τr)^3)+2 (E^((2 tC τc+tCC τc+2 tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (((E0-E∞))^2) (J0-J∞) ν τc ((τr)^3)-(E^((tBB[ec1] τc+tC τc+tCC τc+tD[ec1] τc+2 tBB[ec1] τr+tD[ec1] τr)/(τc τr))) (ec1+ν (tAA-tBB[ec1]+τc)) (((τc-τr))^2) (E∞ τc-E0 τr)-2 (E^((tCC τc+2 tD[ec1] τc+2 tBB[ec1] τr+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)-2 (E^((2 tBB[ec1] τc+tD[ec1] τc+tCC τr+tD[ec1] τr+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)+2 (E^((2 tD[ec1] τc+tCC τr+tBB[ec1] (τc+τr)+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)+2 (E^((tCC τc+tD[ec1] τc+tD[ec1] τr+tBB[ec1] (τc+τr)+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)-(E^((tC τc+tCC τc+2 tD[ec1] τc+tB[ec1] τr+2 tBB[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)-(E^((2 tBB[ec1] τc+tC τc+tD[ec1] τc+tB[ec1] τr+tCC τr+tD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)+(E^((tC τc+2 tD[ec1] τc+tB[ec1] τr+tCC τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)+(E^((tC τc+tCC τc+tD[ec1] τc+tB[ec1] τr+tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)-2 (E^((2 tC τc+tCC τc+tD[ec1] τc+2 tBB[ec1] τr+tD[ec1] τr)/(τc τr))) (E0-E∞) ν ((τr)^2) (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2))+2 (E^((2 tC τc+tD[ec1] τc+tCC τr+tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) ν ((τr)^2) (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2))+(E^((tB[ec1] τc+2 tBB[ec1] τc+tC τc+tCC τr+2 tD[ec1] τr)/(τc τr))) (((E0-E∞))^2) (J0-J∞) τc ((τr)^2) (ec1+ν (-tB[ec1]+τr))-(E^((tB[ec1] τc+tC τc+tCC τc+2 tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (((E0-E∞))^2) (J0-J∞) τc ((τr)^2) (ec1+ν (-tB[ec1]+τr))+(E^((tB[ec1] τc+tC τc+tCC τc+tD[ec1] τc+2 tBB[ec1] τr+tD[ec1] τr)/(τc τr))) (E0-E∞) τr (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2)) (ec1+ν (-tB[ec1]+τr))-(E^((tB[ec1] τc+tC τc+tD[ec1] τc+tCC τr+tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) τr (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2)) (ec1+ν (-tB[ec1]+τr))-(E^((2 tBB[ec1] τc+tC τc+tD[ec1] τc+tBB[ec1] τr+tCC τr+tD[ec1] τr)/(τc τr))) (E0-E∞) (((τc-τr))^2) τr (ec1+ν (tAA-tBB[ec1]+τr))-(E^((tC τc+tCC τc+2 tD[ec1] τc+2 tBB[ec1] τr+tD[ec1] τr)/(τc τr))) (E0-E∞) (τc-τr) τr (ec1 (τc+E∞ J0 τc-E∞ J∞ τc-τr)+ν (-E∞ J0 ((τc)^2)+E∞ J∞ ((τc)^2)-2 tC (τc+E∞ J0 τc-E∞ J∞ τc-τr)+tD[ec1] (τc+E∞ J0 τc-E∞ J∞ τc-τr)-τc τr+E0 J0 τc τr-E∞ J0 τc τr-E0 J∞ τc τr+E∞ J∞ τc τr+((τr)^2)))+(E^((tC τc+2 tD[ec1] τc+tCC τr+tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) (τc-τr) τr (ec1 (τc+E∞ J0 τc-E∞ J∞ τc-τr)+ν (-E∞ J0 ((τc)^2)+E∞ J∞ ((τc)^2)-2 tC (τc+E∞ J0 τc-E∞ J∞ τc-τr)+tD[ec1] (τc+E∞ J0 τc-E∞ J∞ τc-τr)-τc τr+E0 J0 τc τr-E∞ J0 τc τr-E0 J∞ τc τr+E∞ J∞ τc τr+((τr)^2)))+(E^((tC τc+tBB[ec1] (τc+τr)+(tCC+tD[ec1]) (τc+τr))/(τc τr))) (((τc-τr))^3) (ec1 E∞+ν (-E0 τr+E∞ (tAA-tCC+τc+τr)))-(E^((2 tBB[ec1] τc+tC τc+tD[ec1] τc+tCC τr+2 tD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (τc-τr) τr (ec1 E∞+ν (E0 τr-E∞ (2 tC-tD[ec1]+τc+τr)))+(E^((tC τc+tCC τc+tD[ec1] τc+2 tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) τc (τc-τr) τr (ec1 E∞+ν (E0 τr-E∞ (2 tC-tD[ec1]+τc+τr))))

ef2[t_,ec1_]:=(1/((1+ec1) (((τc-τr))^3)))(E^(-((tC τc+tCC τc+tD[ec1] τc+tDD[ec1] τc+t τr+tBB[ec1] τr+tD[ec1] τr)/(τc τr)))) (J0-J∞) (2 (E^((tC τc+2 tCC τc+tD[ec1] τc+tBB[ec1] τr+tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) ν (((τc-τr))^2) ((τr)^2)+2 (E^((tBB[ec1] τc+2 tC τc+tDD[ec1] τc+tCC τr+2 tD[ec1] τr)/(τc τr))) (((E0-E∞))^2) (J0-J∞) ν τc ((τr)^3)-2 (E^((tBB[ec1] τc+2 tC τc+tCC τc+2 tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (((E0-E∞))^2) (J0-J∞) ν τc ((τr)^3)-2 (E^((2 tD[ec1] τc+tDD[ec1] τc+tBB[ec1] τr+tCC τr+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)+2 (E^((tBB[ec1] τc+tD[ec1] τc+tDD[ec1] τc+tCC τr+tD[ec1] τr+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)+2 (E^((tCC τc+2 tD[ec1] τc+tBB[ec1] τr+tDD[ec1] τr+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)-2 (E^((tBB[ec1] τc+tCC τc+tD[ec1] τc+tD[ec1] τr+tDD[ec1] τr+tC (τc+τr))/(τc τr))) (E0-E∞) (J0-J∞) ν ((τc)^2) τr (E∞ τc-E0 τr)-(E^((tC τc+2 tD[ec1] τc+tDD[ec1] τc+tB[ec1] τr+tBB[ec1] τr+tCC τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)+(E^((tBB[ec1] τc+tC τc+tD[ec1] τc+tDD[ec1] τc+tB[ec1] τr+tCC τr+tD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)+(E^((tC τc+tCC τc+2 tD[ec1] τc+tB[ec1] τr+tBB[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)-(E^((tBB[ec1] τc+tC τc+tCC τc+tD[ec1] τc+tB[ec1] τr+tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (ec1+ν (-tB[ec1]+τc)) τr (-E∞ τc+E0 τr)+2 (E^((2 tC τc+tCC τc+tD[ec1] τc+tBB[ec1] τr+tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) ν ((τr)^2) (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2))-2 (E^((2 tC τc+tDD[ec1] τc+tBB[ec1] τr+tCC τr+tD[ec1] (τc+τr))/(τc τr))) (E0-E∞) ν ((τr)^2) (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2))-(E^((tB[ec1] τc+tBB[ec1] τc+tC τc+tDD[ec1] τc+tCC τr+2 tD[ec1] τr)/(τc τr))) (((E0-E∞))^2) (J0-J∞) τc ((τr)^2) (ec1+ν (-tB[ec1]+τr))+(E^((tB[ec1] τc+tBB[ec1] τc+tC τc+tCC τc+2 tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (((E0-E∞))^2) (J0-J∞) τc ((τr)^2) (ec1+ν (-tB[ec1]+τr))+(E^((tB[ec1] τc+tC τc+tD[ec1] τc+tDD[ec1] τc+tBB[ec1] τr+tCC τr+tD[ec1] τr)/(τc τr))) (E0-E∞) τr (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2)) (ec1+ν (-tB[ec1]+τr))-(E^((tB[ec1] τc+tC τc+tCC τc+tD[ec1] τc+tBB[ec1] τr+tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) τr (((τc)^2)+(-2-E∞ J0+E0 (J0-J∞)+E∞ J∞) τc τr+((τr)^2)) (ec1+ν (-tB[ec1]+τr))+(E^((tC τc+tD[ec1] τc+tDD[ec1] τc+tCC τr+tD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) (((τc-τr))^2) τr (ec1+ν (tAA-tBB[ec1]+τr))-(E^((tC τc+tCC τc+tD[ec1] τc+tD[ec1] τr+tDD[ec1] τr+tBB[ec1] (τc+τr))/(τc τr))) (E0-E∞) (((τc-τr))^2) τr (ec1+ν (tAA-tBB[ec1]+τr))+(E^((tC τc+tCC τc+tD[ec1] τc+tDD[ec1] τc+tBB[ec1] τr+tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (((τc-τr))^3) (ec1 E∞+ν (E∞ (tAA-2 tCC+tDD[ec1]-τc-τr)+E0 τr))-(E^((tC τc+2 tD[ec1] τc+tDD[ec1] τc+tBB[ec1] τr+tCC τr+tD[ec1] τr)/(τc τr))) (E0-E∞) (τc-τr) τr (ec1 (τc+E∞ J0 τc-E∞ J∞ τc-τr)+ν (-E∞ J0 ((τc)^2)+E∞ J∞ ((τc)^2)-2 tC (τc+E∞ J0 τc-E∞ J∞ τc-τr)+tD[ec1] (τc+E∞ J0 τc-E∞ J∞ τc-τr)-τc τr+E0 J0 τc τr-E∞ J0 τc τr-E0 J∞ τc τr+E∞ J∞ τc τr+((τr)^2)))+(E^((tC τc+tCC τc+2 tD[ec1] τc+tBB[ec1] τr+tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) (τc-τr) τr (ec1 (τc+E∞ J0 τc-E∞ J∞ τc-τr)+ν (-E∞ J0 ((τc)^2)+E∞ J∞ ((τc)^2)-2 tC (τc+E∞ J0 τc-E∞ J∞ τc-τr)+tD[ec1] (τc+E∞ J0 τc-E∞ J∞ τc-τr)-τc τr+E0 J0 τc τr-E∞ J0 τc τr-E0 J∞ τc τr+E∞ J∞ τc τr+((τr)^2)))+(E^((tBB[ec1] τc+tC τc+tD[ec1] τc+tDD[ec1] τc+tCC τr+2 tD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (τc-τr) τr (ec1 E∞+ν (E0 τr-E∞ (2 tC-tD[ec1]+τc+τr)))-(E^((tBB[ec1] τc+tC τc+tCC τc+tD[ec1] τc+2 tD[ec1] τr+tDD[ec1] τr)/(τc τr))) (E0-E∞) (J0-J∞) τc (τc-τr) τr (ec1 E∞+ν (E0 τr-E∞ (2 tC-tD[ec1]+τc+τr)))+(E^((tC τc+tD[ec1] τc+tDD[ec1] τc+tBB[ec1] τr+tD[ec1] τr+tCC (τc+τr))/(τc τr))) (((τc-τr))^2) (ec1 E∞ (-τc+τr)+ν (-E0 τr (τc+τr)+E∞ (-tAA τc+tCC τc+((τc)^2)+tAA τr-tCC τr+((τr)^2)))))

Fibril strain

ef[t_,ec1_]:=If[t<tB[ec1],0,If[t<tC,1/(1+ec1)(ee1[t]-ec1),If[t<tD[ec1],1/(1+ec1)(ee2[t]-ec1),If[t<tBB[ec1],e11[t,ec1]+e12[t,ec1],If[t<tCC,1/(1+ec1)(ee3[t]-ec1),If[t<tDD[ec1],1/(1+ec1)(ee4[t]-ec1),e11[t,ec1]+e12[t,ec1]+ef1[t,ec1]+ef2[t,ec1]]]]]]]

Plots of imposed strain, fibril strain and fibril stress

Plot[100e[t],{t,0,2},AspectRatio->0.3,PlotStyle->Thick,AxesLabel->{t,"ϵ(%)"}]

Plot[100ef[t,0.03],{t,0,2},AspectRatio->0.3,PlotStyle->Thick,AxesLabel->{t,"ef(%)"},PlotRange->All]

Plot[σf[t,0.03],{t,0,2},AspectRatio->0.3,PlotStyle->Thick,AxesLabel->{t,"σf"}]

Calculate change in fibril length distribution

(*current fibril length divided by initial fascicle length is given by
(l+Δl)/L=(1+ϵc)(1+ϵf)*)

data=RandomVariate[NormalDistribution[7/100,3/100],1000]; (*1000 fibrils sampled from a critical strain distribution with μ=0.07, σ=0.03*)

deformed[t_,i_]:=If[data[[i]]<emax,(1+data[[i]])(1+ef[t,data[[i]]]),1+data[[i]]](*calculate deformed fibril lengths using equation given above*)

SmoothHistogram[Table[deformed[0,i],{i,1,Length[data]}]](*plot undeformed fibril length distribution*)

dist[t_]:=Quiet[SmoothHistogram[{data+1,Table[deformed[t,i],{i,1,Length[data]}]},PlotRange->{0,45},Ticks->{True,False}]](*plot histogram of deformed fibril length distribution at time t*)

distB[t_]:=Show[dist[t],AxesLabel->{"fibril length/fascicle length","frequency density"}] (*as above with axes labelled*)

distsB=Table[distB[t],{t,0,2,0.1}];(*list of length distributions for time between 0s and 2s in increments of 0.1s*)

time=Table[Plot[{e[t],1000(t- T)},{t,0,2},PlotRange->{{-0.01,2.01},{0,0.06}},AxesLabel->{"time","strain"}],{T,0,2,0.1}];(*list of images showing current time against imposed strain*)

ListAnimate[Table[GraphicsColumn[{time[[Ceiling[i/5]]],distsB[[Ceiling[i/5]]]}],{i,1,5Length[time]}],1](*animate previously calculated images*)

distsC=Table[GraphicsColumn[{time[[Ceiling[i/5]]],distsB[[Ceiling[i/5]]]}],{i,1,5Length[time]}];

Export["INSERT FILEPATH",distsC]

What happens if we neglect fibril creep and impose fibril strain for all time after the initial taughtening?

efNoCreep[t_,ec1_]:=If[t<tB[ec1],0,1/(1+ec1)(e[t]-ec1)](*the fibril strain that results from the imposed strain if creep is nelgected*)

deformedNoCreep[t_,i_]:=If[data[[i]]<emax,(1+data[[i]])(1+efNoCreep[t,data[[i]]]),1+data[[i]]](*calculate deformed fibril lengths using equation given above*)

distNoCreep[t_]:=Quiet[SmoothHistogram[{data+1,Table[deformedNoCreep[t,i],{i,1,Length[data]}]},PlotRange->{0,45},Ticks->{True,False}]](*plot histogram of deformed fibril length distribution at time t*)

distBNoCreep[t_]:=Show[distNoCreep[t],AxesLabel->{"fibril length/fascicle length","frequency density"}] (*as above with axes labelled*)

distsBNoCreep=Table[distBNoCreep[t],{t,0,2,0.1}];(*list of length distributions for time between 0s and 2s in increments of 0.1s*)

ListAnimate[Table[GraphicsColumn[{time[[Ceiling[i/5]]],distsBNoCreep[[Ceiling[i/5]]]}],{i,1,5Length[time]}],1](*animate previously calculated images*)

distsCNoCreep=Table[GraphicsColumn[{time[[Ceiling[i/5]]],distsBNoCreep[[Ceiling[i/5]]]}],{i,1,5Length[time]}];

Export["INSERT FILEPATH",distsCNoCreep]