// 1 declare the main function func sphere (sigma=, DCa=, buffers=, rmin=, rmax=, xlog=, ylog=, nr=, nstep=, t0=, trun=, bc0=, bc1=, prompt=, tol=) { // 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 (rmin)) rmin = 0.01; if (is_void (rmax)) rmax = 2.0; if (is_void (nr)) nr = 30; if (is_void (nstep)) nstep = 100; if (is_void (trun)) trun = 10; if (is_void (tol)) tol = 1.e-8; if (is_void (bc0)) bc0 = 0; if (is_void (bc1)) bc1 = 0; if (is_void (xlog)) xlog = 1; if (is_void (ylog)) ylog = 1; if (is_void (prompt)) prompt = 0; if (is_void (t0)) t0 = trun * 1.e-5; // 3 for convenience, if there are no buffers, add a dummy one // which doesnt bind calcium if (is_void (buffers)) buffers = ([[1., 0, 0, 1]]); // 4 set the time points at which to evaluate the solution tpts = grow (0., exp (span(log(t0), log(trun), nstep))); // 5 divide the solution domain into nr logarithmic shells rpts = grow (0., exp (span (log(rmin), log(rmax), nr))); // 6 find the centers of the volumes between shells; rcen = rpts(zcen); // 7 find the volumes of the compartments vol = (4./3. * pi * rpts^3)(dif); // 8 the diffusion rate between compartments is given by the area // connecting them divided by the distance between them rdif = (4 * pi * rpts^2) (2:nr) / (rpts(zcen))(dif); // 9 set up the initial conditions. np = numberof (vol); nbuf = numberof (buffers(1,)); conc = array (double, 1 + 2 * nbuf, np); conc(1,) = 0.; conc (2*indgen(nbuf),) = buffers(1,); // 10 if bc0 is 1, then the concnetration at the centere is fixed at // sigma, instead of injecting calcium at rate sigma. if (bc0 == 1) conc(1,1) = sigma; // 11 convert the set of buffers into a general reaction matrix; // rreac will contain the reaction rates, ireac a list of // reactants and products for each reaction, and mobil the // mobility of each species rreac = []; ireac = []; mobil = [DCa]; for (i = 1; i <= nbuf; i++) { // the forward reaction Ca + B -> BCa grow, rreac, buffers(2, i); grow, ireac, [[1, 2*i, 0, 2*i+1, 0, 0]]; // reverse reaction BCa -> B + Ca grow, rreac, buffers(3, i); grow, ireac, [[2*i+1, 0, 0, 1, 2*i, 0]]; grow, mobil, [buffers(4,i), buffers(4,i)]; } // 12 now run the system for nstep timesteps; for (istep = 1; istep <= nstep; istep++) { dt = tpts(istep+1) - tpts(istep); iit = 0; // 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, rdif, vol, ireac, rreac, dt, sigma, bc0, bc1); // 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, rdif, vol, ireac, rreac, dt, bc0, bc1); // 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, rdif, vol, ireac, rreac, dt, sigma, bc0, bc1); // 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; // 21 draw a graph showing the various concentrations fma; rp = (xlog == 1 ? log10(rcen) : rcen); cp = (ylog == 1 ? log10(conc + 1.e-12) : conc); plg, cp(1,), rp, marks=0, color = "white"; for (i = 1; i <= nbuf; i++) { plg, cp (2*i,), rp, marks=0, color = "blue"; plg, cp (2*i+1,), rp, marks=0, color = "black", type=2; } if (xlog == 1) { xytitles, "log10 r", "log10 concentration"; } else { xytitles, "r", "concentration"; } pltitle, swrite (format = "concentration at time %8.5f", tpts(istep+1)); // 22 if the prompt keyword is set, wait for the user to // press return before continuing if (prompt) read, aaa; } } func rdEquns (conc, dconc, mobil, rdif, vol, ireac, rreac, dt, sigma, bc0, bc1) { 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(4) > 0) R(ir(4),) -= r; if (ir(5) > 0) R(ir(5),) -= 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 * (rdif / vol(1:np-1))(-,); R(,2:np) += deltc * mobil(,-) * dt * (rdif / vol(2:np))(-,); // 34 the bc0 and bc1 flags set what type of boundary conditions // 0 means free, 1 fixed // for the center, free means sigma is taken as an injection rate // fixed means the concentration stays fixed at sigma if (bc0 == 0) { R(1,1) -= sigma * dt / vol(1); } else { R(1,1) = dconc(1,1); } // 35 for the outer boundary, fixed means the concnetration is kept at its // start value as though there were an infinite sink or source if (bc1 == 1) R(,np) = dconc(,np); return R; } func rdDerivs (conc, dconc, mobil, rdif, vol, ireac, rreac, dt, bc0, bc1) { /* 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(4); if (m > 0) M(m,k,,2) -= r; m = ir(5); 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(4); if (m > 0) { M(m, k1,,2) -= r * conc(k2,); M(m, k2,,2) -= r * conc(k1,); } m = ir(5); if (m > 0) { M(m, k1,,2) -= r * conc(k2,); M(m, k2,,2) -= r * conc(k1,); } } } // derivatives of diffusion terms; fa = dt * mobil(,-) * (rdif / vol(1:np-1))(-,); fb = dt * mobil(,-) * (rdif / vol(2:np))(-,); 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,); } if (bc0 == 1) { M(1,,1,2:3) = 0.; M(1,1,1,2) = -1.; } if (bc1 == 1) { M(,,np,1:2) = 0.; for (i = 1; i <= numberof (M(,1,1,1)); i++) M(i,i,np,2) = -1.; } 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; } func tstMat(dum) { ne = 6; nb = 7; M = random (ne, ne, nb, 3); R = random (ne, nb); rf = R(*); mf = array (double, ne*nb, ne*nb); for (i = 1; i <= nb; i++) { a = ne * (i-1); if (i > 1) mf(a+1:a+ne, a-ne+1:a) = M(,,i,1); mf(a+1:a+ne, a+1:a+ne) = M(,,i,2); if (i < nb) mf(a+1:a+ne, a+ne+1:a+ne+ne) = M(,,i,3); } z = LUsolve (mf, rf); W = tbdMatrixSoln (M, R); max (abs (W(*) - z)); }