Matlab code 4: Spiro's model for chemotaxis Filename: spiro.m clear; close; To=8e-6; Yo=20e-6; Bo=1.7e-6; options = odeset('RelTol',1e-9,'AbsTol',[1e-9 1e-9 1e-9 1e-9 1e-9]); [t y]=ode23('spirofunc',[0 80],[4e-6 4e-6 0e-6 0.5e-6 10e-6],options); Ptot=To-y(:,1)-y(:,2); Bp=y(:,4); Yp=y(:,5); metlevel=1-(y(:,1)+y(:,3))/To; phoslevel=1-(y(:,1)+y(:,2))/To; subplot(2,2,1) plot(t,phoslevel,'bx'); axis([10 80 0 0.1]); title('Phosphorylation Level'); subplot(2,2,2) plot(t,metlevel,'rx'); title('Methylation Level'); axis([10 80 0 1]); subplot(2,2,3) plot(t,Bp/Bo,'gx'); axis([10 80 0 1]); title('Bp/Btot'); subplot(2,2,4) plot(t,Yp/Yo,'yx'); axis([10 80 0 1]); title('Yp/Ytot'); Filename: spirofunc.m function dydt = f(t,y,flag) % constants from Table 3 (Spiro et al.) k1c=0.17; % 1/s k3c=30*k1c; % 1/s ratiok1bk1a=1.7e-6; % M ratiok3ck3a=1.7e-6; % M k_1=4e5; % 1/(Ms) k_3=k_1; % 1/(Ms)| k8=15; % 1/s k9=3*k8; % 1/s k11=0; % 1/s k12=1.1*k8; % 1/s kb=8e5; % 1/(Ms) ky=3e7; % 1/(Ms) k_b=0.35; % 1/s k_y=5e5; % 1/(Ms) Kbind=1e6; % 1/M Yo=20e-6; % M Bo=1.7e-6; % M To=8e-6; % M Ro=0.3e-6; % M Zo=40e-6; % M % [T2]+[LT2] = y(1) % [T3]+[LT3] = y(2) % [T2p]+[LT2p] = y(3) % [Bp] = y(4) % [Yp] = y(5) cligand=1e-6; if t>20 cligand=1e-3; end; if t>50 cligand=1e-6; end; Vmaxunbound=k1c*Ro; % maximum turnover rate (MM kinetics) for unbound receptors Vmaxbound=k3c*Ro; % maximum turnover rate (MM kinetics) for bound receptors KR=ratiok1bk1a; % Michaelis constant fb=Kbind*cligand/(1+Kbind*cligand); % fraction receptors bound to ligand fu=1-fb; % fraction receptors not bound to ligand kpt=ky*(Yo-y(5))+kb*(Bo-y(4)); ydot1=(-k8*fu-k11*fb)*y(1)+kpt*y(3)+(k_1*fu+k_3*fb)*y(2)*y(4)-Vmaxunbound*y(1)*fu/(KR+y(1)*fu)-Vmaxbound*y(1)*fb/(KR+y(1)*fb); ydot2=(-k9*fu-k12*fb)*y(2)+kpt*(To-y(1)-y(2)-y(3))-(k_1*fu+k_3*fb)*y(2)*y(4)+Vmaxunbound*y(1)*fu/(KR+y(1)*fu)+Vmaxbound*y(1)*fb/(KR+y(1)*fb); ydot3=(k8*fu+k11*fb)*y(1)-kpt*y(3)+(k_1*fu+k_3*fb)*(To-y(1)-y(2)-y(3))*y(4)-Vmaxunbound*y(3)*fu/(KR+y(3)*fu)+Vmaxbound*y(3)*fb/(KR+y(3)*fb); ydot4=kb*(To-y(1)-y(2))*(Bo-y(4))-k_b*y(4); ydot5=ky*(To-y(1)-y(2))*(Yo-y(5))-k_y*y(5)*Zo; dydt=[ydot1; ydot2; ydot3; ydot4; ydot5];