func tspark (dum) { spark, c0 = [0., 40, 1., 0., 0.], ireac = [[1,0,2,0], [3,1,4,0], [4,0,3,1], [1,4,2,4], [2,4,1,4], [4,0, 5, 0]], rreac = [0.1, 1., 0.1, 0.3, 0.3, 1.1]; } func tfdf (dum) { spark, c0=[0.,40], ireac=[1,0,2,0], rreac=[0.1], fdf=0.01; } // 1 declare the main function func spark (DCa=, sigma=, buffers=, dx=, nsite=, npps=, nstep=, trun=, prompt=, logy=, tol=, ireac=, rreac=, c0=, bufa=, bufb=, nr=, rmax=, fdf=) { // 2 default values for things which haven't been specified in the command if (is_void (sigma)) sigma = 5.; if (is_void (DCa)) DCa = 0.25; if (is_void (rmax)) rmax = 20.0; if (is_void (nstep)) nstep = 400; if (is_void (trun)) trun = 40; if (is_void (tol)) tol = 1.e-8; if (is_void (logy)) logy = 0; if (is_void (prompt)) prompt = 0; if (is_void (ireac)) ireac=[1, 0, 2, 0]; if (is_void (rreac)) rreac=[0.1]; if (is_void (c0)) c0 = [0., 100.]; if (is_void (fdf)) fdf = 0.; if (!is_void (nr)) { if (is_void (rmax)) rmax = 20; } else { if (is_void (nsite)) nsite = 6; if (is_void (dx)) dx = 0.5; if (is_void (npps)) npps = 6; } dt = double(trun) / nstep; // 4 an array of boundaries of cells; if (is_void(nr)) { rpts = span (0., dx*nsite, 1 + npps * nsite); } else { rpts = span (0., rmax, nr+1); } // 5 irs are the indices of the release sites: the points which // will have non-zero channel density; if (is_void(nr)) irs = npps/2 + npps * (indgen(nsite)-1); // 6 find the centers of cells; rcen = rpts(zcen); // 7 count buffers etc nc = numberof (c0); na = (is_array(bufa) ? numberof(bufa(1,)) : 0); nb = (is_array(bufb) ? numberof(bufb(1,)) : 0); ns = nc + 2 * na + 2 * nb; np = numberof (rcen); // 9 set up the initial conditions. conc = array (double, nc+2*na+2*nb, np); conc(1:nc,) = c0; if (na > 0) conc (nc + 2*indgen(na),) = bufa(1,); if (nb > 0) conc (nc+2*na + 2*indgen(nb),) = bufb(1,); // 10 if a channel density has been specified, only apply it at the // release sites, making the rest 0 if (is_void (nr) && numberof(c0) >= 3) { conc(3,) = 0.; conc(3,irs) = c0(3); } conc(1,1) = sigma; // 11 append the buffering reactions to the reaction scheme; mobil = 0. * c0; mobil(0) = mobil(1) = DCa; for (i = 1; i <= na; i++) { grow, rreac, bufa(2, i); grow, ireac, [[nc+1, nc+2*i, nc+2*i+1, 0]]; grow, rreac, bufa(3, i); grow, ireac, [[nc+2*i+1, 0, nc+1, nc+2*i]]; grow, mobil, [buffers(4,i), buffers(4,i)]; } for (i = 1; i <= nb; i++) { grow, rreac, bufb(2, i); grow, ireac, [[nc+2*na+2, nc+2*na+2*i, nc+2*na+2*i+2, 0]]; grow, rreac, bufa(3, i); grow, ireac, [[nc+2*a+2*i+2, 0, nc+2*a+2, nc+2*a+2*i]]; grow, mobil, [buffers(4,i), buffers(4,i)]; } // 11a an empty array to accumulate all the calcium concentrations; cca = []; window, 1; animate, 1; // 11b an array to record which sites have already fired; fireable = array (0, numberof (conc(1,))); fireable(irs) = 1; // 12 now run the system for nstep timesteps; time = 0.0; for (istep = 1; istep <= nstep; istep++) { iit = 0; time += dt; // 13 dconc is the current guess at the changes in concentration // over one step dconc = 0. * conc; // 14 evaluate the equations for this step R = rdEquns (conc, dconc, mobil, ireac, rreac, dt); // 15 use a Newton-Raphson method to iterate towards a solution for // the concentration changes over this step, continuing // until the average of R is below the allowed tolerance for (iit = 0; iit < 12 & avg(abs(R)) > tol; iit++) { // 16 evaluate the matrix of derivatives B = rdDerivs (conc, dconc, mobil, ireac, rreac, dt); // 17 solve a matrix eqution for hte corrections W and apply them to // dconc W = tbdMatrixSoln (B, R); dconc -= W; // 18 evaluate the equations at the new guess, and repeat as // necessary R = rdEquns (conc, dconc, mobil, ireac, rreac, dt); // 19 if its taking too many iterations, something must have gone // wrong if (iit > 6) print, istep, iit, avg(abs(R)); } // 20 increment the concentrations for this timestep conc += dconc; // 20a if doinga fire-diffuse-fire simulation, check whether // any release sites have fired and if so dump their calcium if (fdf > 0. & anyof (conc(1,) * fireable > fdf)) { w = where (conc(1,) * fireable > fdf); conc(1:2,w) = ( (conc(1:2,w))(avg,) )(-,); fireable(w) = 0; } // 21 draw a graph showing the various concentrations window, 0; fma; cp = (logy == 1 ? log10(conc + 1.e-12) : conc); plg, cp(1,), rcen, marks=0, color = "white"; plg, cp(2,), rcen, marks=0, color = "white", type = 2; for (i = 3; i <= nc; i++) { plg, cp(i,), rcen, marks=0, color = "red", type = i-2; } for (i = 1; i <= na; i++) { plg, cp (nc + 2*i,), rp, marks=0, color = "blue"; plg, cp (nc + 2*i+1,), rp, marks=0, color = "blue", type=2; } for (i = 1; i <= nb; i++) { plg, cp (nc + 2*na + 2*i,), rp, marks=0, color = "black"; plg, cp (nc + 2*na + 2*i+1,), rp, marks=0, color = "black", type=2; } xytitles, "r", (logy==1 ? "log10 concentration" : "concentration"); pltitle, swrite (format = "concentration at time %8.5f", time); // 21a accumulate the calcium concentration grow, cca, [conc(1,)]; window, 1; fma; pli, cca, 0, 0, rpts(0), time; // 22 if the prompt keyword is set, wait for the user to // press return before continuing if (prompt) read, aaa; } window, 1; animate, 0; } func rdEquns (conc, dconc, mobil, ireac, rreac, dt) { nreac = numberof(ireac(1,)); // 31 first term in an eqution is just dconc R = dconc; // 32 the reactions terms have the form // R = change in concentration - sum of inputs + sum of outputs // when R = 0 we have a solution ofr the timestep ct = conc + dconc; for (i = 1; i <= nreac; i++) { r = dt * rreac(i) + 0. * (conc(1,)); ir = ireac (,i); if (ir(1) > 0) r *= ct (ir(1),); if (ir(2) > 0) r *= ct (ir(2),); if (ir(1) > 0) R(ir(1),) += r; if (ir(2) > 0) R(ir(2),) += r; if (ir(3) > 0) R(ir(3),) -= r; if (ir(4) > 0) R(ir(4),) -= r; } // 33 the diffusion terms are similar, except here the source // terms are by diffusion from neighbouring points np = numberof (conc(1,)); deltc = ct(,dif); R(,1:np-1) -= deltc * mobil(,-) * dt; R(,2:np) += deltc * mobil(,-) * dt; return R; } func rdDerivs (conc, dconc, mobil, ireac, rreac, dt) { /* this routine fills tri block diagonal a matrix of derivatives of the equations above. The structure of the matrix is in blocks as M(,,1,2) M(,,1,3) 0 0 0 M(,,2,1) M(,,2,2) M(,,2,3) 0 0 0 M(,,3,1) M(,,3,2) M(,,3,3) 0 0 0 M(,,4,1) M(,,4,2) M(,,4,3) 0 0 0 M(,,5,1) M(,,5,2) */ ns = numberof (conc(,1)); np = numberof (conc(1,)); M = array (double, ns, ns, np, 3); for (i = 1; i <= ns; i++) M(i,i,,2) = 1.; // reaction terms; nreac = numberof (rreac); for (i = 1; i <= nreac; i++) { r = dt * rreac(i) + 0. * (conc(1,)); ir = ireac (,i); if (ir(2) == 0) { // spontaneous decay; k = ir(1); M(k, k,,2) += r; m = ir(3); if (m > 0) M(m,k,,2) -= r; m = ir(4); if (m > 0) M(m,k,,2) -= r; } else { // binding reactions; k1 = ir(1); k2 = ir(2); M(k1, k1,,2) += r * conc(k2,); M(k2, k2,,2) += r * conc(k1,); M(k1, k2,,2) += r * conc(k1,); M(k2, k1,,2) += r * conc(k2,); m = ir(3); if (m > 0) { M(m, k1,,2) -= r * conc(k2,); M(m, k2,,2) -= r * conc(k1,); } m = ir(4); if (m > 0) { M(m, k1,,2) -= r * conc(k2,); M(m, k2,,2) -= r * conc(k1,); } } } // derivatives of diffusion terms; fa = dt * mobil(,-); fb = dt * mobil(,-); for (i = 1; i <= ns; i++) { M(i, i, 1:np-1, 2) += fa(i,); M(i, i, 1:np-1, 3) -= fa(i,); M(i, i, 2:np, 2) += fb(i,); M(i, i, 2:np, 1) -= fb(i,); } return M; } func tbdMatrixSoln (M, R) { // 51 solve the matrix equation as described in rdDerivs // the contraction operations make yorick particularly good for // this sort of thing A = M(,,,1); B = M(,,,2); C = M(,,,3); nb = numberof (A(1,1,)); // 52 eliminate the left hand block, diagonalise the central one for (k = 1; k <= nb; k++) { if (k > 1) { B(,,k) -= A(,+,k) * C(+,,k-1); R(,k) -= A(,+,k) * R(+,k-1); } s = LUsolve (B(,,k)); R(,k) = s(,+) * R(+,k); C(,,k) = s(,+) * C(+,,k); } // 53 backsubstitution W = 0. * R; W(,nb) = R(,nb); for (k = nb-1; k >= 1; k--) W(,k) = R(,k) - C(,+,k) * W(+,k+1); return W; }