clear; close; n=100; % number of discretized spatial points L=1; % spatial length dx=L/n; % incremental change in x coordinate dt=0.05; % incremental change in time T=1000; % total integration time lam0=0.1; % constant background level of ligand lamA=0.9; % amplitude of ligand peak tplot=1; % plot frequency A=3; % magnitude of peak at trailing edge (choose 0 for uniform distribution) %rom Table 1 (Narang et al.) kappaeps=0.02; kappa_f=1.1; kappa_i=0.1; kappa_p=0.1; eps_d=0.1; psi_p=0.01; psi_i=0.01; delta_p=0.0001; delta_r=0.0001; % initial conditions pold=0.0991*ones(1,n); iold=0.0991; sigma=0.5*n; for x=1:n, lam(x)=lam0*(1+lamA*exp(-(x-n/2)^2/sigma)); lamold(x)=lam0*(1+lamA*exp(-(x-n/2)^2/sigma)); end; rhoold=(1+1/lam0+eps_d)./(1+1./lam+eps_d); % plot initial conditions xvec=linspace(0,1,n); subplot(2,2,1); plot(xvec,lam,'bo'); title('cAMP GRADIENT'); subplot(2,2,2); plot(xvec,rhoold,'go'); axis([0 1 0 2]); title('ACTIVE RECEPTORS'); subplot(2,2,3); plot(xvec,pold,'ro'); title('PIP2s (ACTIVATOR)'); axis([0 1 0 1]); subplot(2,2,4); ivec=iold*ones(n,1); plot(xvec,ivec,'yo'); title('IP3 (INHIBITOR)'); axis([0 1 0 0.2]); drawnow; input('Hit a key to start simulation'); close; figure; % start of iterative time loop for timestep=1:(T/dt); t=timestep*dt; if t==50 input('Hit a key to apply ligand distribution at trailing edge'); close; figure; for x=1:n/2, lam(x)=lam0*(1+A*exp(-(x)^2/sigma)); lamnew(x)=lam(x); end; for x=n/2+1:n, lam(x)=lam0*(1+A*exp(-(x-n)^2/sigma)); lamnew(x)=lam(x); end; rhonew=rhoold.*(1+1./lamold+eps_d)./(1+1./lamnew+eps_d); rhoold=rhonew; end; % equations for x=1, incorporating periodic boundary conditions drhodt(1)=(1-rhoold(1))*kappaeps/(1+1/lam(1)+eps_d) + ... delta_r*(rhoold(n)-2*rhoold(1)+rhoold(2))/dx^2; dpdt(1)=kappa_f*rhoold(1)*pold(1)^2*(1-pold(1))-pold(1)*iold + ... psi_p-kappa_p*pold(1)+ delta_p*(pold(n)-2*pold(1)+pold(2))/dx^2; % equation for internal points (from 2 to n-1) for x=2:n-1, drhodt(x)=(1-rhoold(x))*kappaeps/(1+1/lam(x)+eps_d) + ... delta_r*(rhoold(x-1)-2*rhoold(x)+rhoold(x+1))/dx^2; dpdt(x)=kappa_f*rhoold(x)*pold(x)^2*(1-pold(x))-pold(x)*iold + ... psi_p-kappa_p*pold(x)+ delta_p*(pold(x-1)-2*pold(x)+pold(x+1))/dx^2; end; % equations for x=n, incorporating periodic boundary conditions drhodt(n)=(1-rhoold(1))*kappaeps/(1+1/lam(1)+eps_d) + ... delta_r*(rhoold(n-1)-2*rhoold(n)+rhoold(1))/dx^2; dpdt(n)=kappa_f*rhoold(n)*pold(n)^2*(1-pold(n))-pold(n)*iold + ... psi_p-kappa_p*pold(n-1)+ delta_p*(pold(n-1)-2*pold(n)+pold(1))/dx^2; % calculate di/dt didt=kappa_f*mean(rhoold)*mean(pold)^2*mean(1-pold)-mean(pold)*iold + psi_i-kappa_i*iold; % update the new values of rho, p, and i. for x=1:n, rhonew(x,timestep)=rhoold(x)+drhodt(x)*dt; pnew(x,timestep)=pold(x)+dpdt(x)*dt; rhoold(x)=rhoold(x)+drhodt(x)*dt; pold(x)=pold(x)+dpdt(x)*dt; end; inew(timestep)=iold+didt*dt; iold=iold+didt*dt; if rem(t,tplot)==0 xvec=linspace(0,1,n); subplot(2,2,1); plot(xvec,lam,'bo'); title('cAMP GRADIENT'); subplot(2,2,2); plot(xvec,rhoold,'go'); title('ACTIVE RECEPTORS'); axis([0 1 0 2]); subplot(2,2,3); plot(xvec,pold,'ro'); title('PIP2s (ACTIVATOR)'); axis([0 1 0 1]); subplot(2,2,4); ivec=iold*ones(n,1); plot(xvec,ivec,'yo'); title('IP3 (INHIBITOR)'); axis([0 1 0 0.2]); drawnow; t iold end; end;