2f:Tf02,{ // generator plotting const time = solution.t; const sol = solution.sol; // Parameters for the Model const colArea = typeParams.colArea; const caDark = typeParams.concCaDark; const cgDark = typeParams.concCgDark; const betaDark = typeParams.betaDark; const betaSub = typeParams.betaSub; const muRa = typeParams.muRa; const muTa = typeParams.muTa; const muPa = typeParams.muPa; const xi = typeParams.xi; // Cyclase parameters const nAlpha = typeParams.nAlpha; const rAlpha = typeParams.rAlpha; const KAlpha = typeParams.KAlpha / caDark; // Exchanger and Channel properties const fChCa = typeParams.fChCa; const nCh = typeParams.nCh; const KEx = typeParams.KEx / caDark; const KCh = typeParams.KCh / cgDark; // Rate for the Calcium dynamics const muCa = typeParams.muCa; const variables = ['Ra', 'Ta', 'Pa', 'y', 'z']; const ind_Ra = variables.indexOf('Ra'); const ind_Ta = variables.indexOf('Ta'); const ind_Pa = variables.indexOf('Pa'); const ind_y = variables.indexOf('y'); const ind_z = variables.indexOf('z'); function getResponse(c,a){ // % cGMP and Calcium normalized by the dark values let cgmp = c.map(v => Math.exp(-v)); let ca = a.map(v => Math.exp(-v)); // % Currents normalized by the dark values // % intracellular recording (whole-cell patch) // % output.curr = 1-output.current; let currCh = cgmp.map(val => 2/(fChCa+2)*(1-(1+KCh**nCh)/(val**nCh + KCh**nCh)*val**nCh)); let currEx = ca.map(val => fChCa / (fChCa + 2) * (1 - (1 + KEx) / (val + KEx) * val)); //let photocurrent = currCh + currEx; return currEx.map((v,i) => (v+currCh[i]) * maxResponse); } // % cGMP and Calcium normalized by the dark values let cgmp = sol[ind_y].map(v => Math.exp(-v)); let ca = sol[ind_z].map(v => Math.exp(-v)); // % Currents normalized by the dark values // % intracellular recording (whole-cell patch) // % output.curr = 1-output.current; let currCh = cgmp.map(val => 2/(fChCa+2)*(1-(1+KCh**nCh)/(val**nCh + KCh**nCh)*val**nCh)); let currEx = ca.map(val => fChCa / (fChCa + 2) * (1 - (1 + KEx) / (val + KEx) * val)); //let photocurrent = currCh + currEx; //let photocurrent = currEx.map((v,i) => (v+currCh[i]) * maxResponse); let photocurrent = getResponse(sol[ind_y],sol[ind_z]); // maxResponse is in pA let d = time[0].map((t,i) => {return {time:t,response:photocurrent[i]}}); // generate a short segment at a time const nChunks = 120; const cSize = time[0].length / nChunks; const tChunks = Array(nChunks) .fill(Array(cSize).fill(0).map((v,i) => i)) .map((v,i) => v.map(t => time[0][t+(i*cSize)])); var sIdx = 0; var seg = tChunks[sIdx]; var cSol = await getChunk(solution,[seg[0],seg.slice(-1)[0]]); var resp = getResponse(cSol.sol[ind_y],cSol.sol[ind_z]); var chunkData = cSol.t[0].map((t,i) => {return {time:t,response:resp[i]}}); var iT = Date.now(); var eT = Date.now(); const dT = cSize/sampleRate * 1000; //ms while (true) { eT = Date.now(); if (eT - iT > dT){ iT = Date.now(); sIdx = (sIdx+1) % nChunks; seg = tChunks[sIdx]; cSol = await getChunk(solution,[seg[0],seg.slice(-1)[0]]); resp = getResponse(cSol.sol[ind_y],cSol.sol[ind_z]); chunkData = cSol.t[0].map((t,i) => {return {time:t,response:resp[i]}}); } yield Plot.plot({ x: { label: "Time (s)" }, y: { label: "Photocurrent (pA)", grid: true }, marginTop: 30, marginBottom: 45, marks: [ Plot.ruleY([0]), Plot.ruleX([0]), Plot.lineY(d, {x: "time", y: "response"}), Plot.lineY(chunkData,{x: "time", y: "response",stroke:'red',strokeWidth:10}) ] }); } }30:T7b4,typeParams = { let params; switch(cellType){ case "Rod": params = { betaDark: 4.1, // 1/s currDark: 15.0, // pA concCaDark: 0.3, // muM concCgDark: 4.0, // muM fChCa: 0.12, colArea: 0.28, betaSub: 0.021, // Parameters for PDE activation xi: 0.45/2.5, muRa: 28., muTa: 23., muPa: 5., // cGMP synthesis nAlpha: 2., KAlpha: 0.14, // muM rAlpha: 1/30, // CNG channel nCh: 2.5, KCh: 20., // muM // Exchanger KEx: 1.6, // muM muCa: 50. // Effective rate for the calcium dynamics }; break; case "Cone": params = { betaDark: 11, // 1/s currDark: 15, // pA concCaDark: 0.3, // muM concCgDark: 4, // muM fChCa: 0.3, colArea: 0.013, betaSub: 0.021, // Parameters for PDE activation xi: 0.0018/2.5, muRa: 70, muTa: 71, muPa: 37.8, // cGMP synthesis nAlpha: 2, KAlpha: 0.14, rAlpha: 1/30, // CNG channel nCh: 2.5, KCh: 20, // muM // Exchanger KEx: 1.6, // muM // Effective rate for the calcium dynamics muCa: 150 } break; } // let Rstar = 1; // # isomerisations // let flashDuration = 0.010; // let fs = 1000; // let responseDuration = 2.0; // let colArea = 0.28; // let flashIntensities = Rstar / colArea / flashDuration; // let flashTimes = [0, flashDuration]; // let time = Array.from({length: responseDuration * fs}, (_, i) => i / fs); let fs = sampleRate; params['Rstar'] = Rstar; params['flashTimes'] = [0,stimulusDuration]; params['flashIntensity'] = Rstar / params.colArea / stimulusDuration; params['time'] = Array.from({length: responseDuration * fs}, (_, i) => i / fs - stimulusOffset); return params }31:T98f,function ntrp15s(tinterp, tnew, ynew, h, dif, k, idxNonNegative) { // NTRP15S Interpolation helper function for ODE15S. // [yinterp, ypinterp] = ntrp15s(tinterp, t, y, tnew, ynew, h, dif, k, idxNonNegative) uses data computed in // ODE15S to approximate the solution at time tinterp. tinterp may be a // scalar or a row vector. // The arguments t and y do not affect the computations. They are required // for consistency of syntax with other interpolation functions. Any values // entered for t and y are ignored. // // [yinterp, ypinterp] = ntrp15s(tinterp, t, y, tnew, ynew, h, dif, k, idxNonNegative) returns // also the derivative of the polynomial approximating the solution. // // idxNonNegative has indices of solution components that must be non-negative. Negative // yinterp(idxNonNegative) are replaced with zeros and the derivative ypinterp(idxNonNegative) is // set to zero. // Requires: math.js@11.8+ const s = (tinterp - tnew)/h; var ypinterp; var yinterp; // becomes math.js matrix const neq = ynew.length; // n rows const nt = Array.isArray(tinterp) ? tinterp.length : 1; // n time points to interp // console.log("ntrp15s"); // console.log('tinterp:',tinterp); // console.log("tnew:",tnew); // console.log("ynew:",ynew); // console.log("h:",h); // console.log("dif:",dif); // console.log("k:",k); // console.log("idxNonNegative:",idxNonNegative); if (k === 0) { yinterp = math.add(ynew, math.multiply(math.column(dif,0), s)); } else { // cumprod collapses vectors const K = math.range(0, k, true).resize([k+1,1]).valueOf(); // column const kI = K.map((v,i) => Array(nt).fill(v[0])); yinterp = math.add( math.subset(ynew,math.index(math.range(0,neq),math.zeros(nt).valueOf())), math.multiply( math.subset(dif,math.index(math.range(0,neq),K.map(k => k[0]))), cumprod(math.chain(math.add(s,kI)).dotDivide(math.add(kI,1)).done()) ) ); } // Non-negative solution if (idxNonNegative.length > 0) { // locate yinterp rows that shouldn't be negative let idx = yinterp.map( (v,i) => idxNonNegative.includes(i[0]) ? i : [false,false] ) .valueOf() .filter(v => v[0][0]) .map(v => v[0]); if (idx.length > 0) { math.apply(idx,1,ii => yinterp.subset(math.index(...ii),math.max(0,yinterp.subset(math.index(...ii))))); } } return yinterp; }32:T64e,function odenumjac(F, t, y, f0, thresh, fac) { // Initialize. const ny = y.length; // y is column array const nF = f0.length; // Fvalue is column matrix if (fac === undefined) { fac = math.matrix(Array(ny).fill(Math.sqrt(Number.EPSILON)).map(v=>[v])); } // Select an increment del for a difference approximation to // column j of dFdy. The vector fac accounts for experience // gained in previous calls to numjac. const yscale = math.resize(math.max(math.concat(math.abs(y),thresh),1),[ny,1]); let del = math.subtract(math.add(y,math.dotMultiply(fac,yscale)),y); if (nF === ny) { const s = math.largerEq(f0,0); del = math.dotMultiply(math.subtract(s,math.not(s)),math.abs(del)); } // Form a difference approximation to all columns of dFdy. const ydel = math.add(y, math.diag(math.squeeze(del))); // use math.js for matrix operations const Fdel = math.zeros(nF, ny); // for each equation, get partial derivative by mapping to ydel let ycol, fdelTmp; for (let j = 0; j < ny; j++) { // F(t0,yest,args) ycol = math.column(ydel,j); fdelTmp = F(t,ycol); // input vector Fdel.subset(math.index(math.range(0,ny),j),math.matrix(fdelTmp).resize([ny,1])); } const nfcalls = ny; // stats const nfevals = ny; // stats (at least one per loop) const Fdiff = math.subtract(Fdel, math.subset(f0,math.index(math.range(0,ny),math.zeros(nF)))); // use math.js for matrix operations const dFdy = math.multiply(Fdiff, math.divide(1, math.diag(math.squeeze(del)))); // use math.js for matrix operations return [dFdy.valueOf(), nfevals]; }33:T4108,async function solve(ode, tspan, y0, options) { // This is a direct port of MATLAB's ode15s function // it isn't terribly fast, unfortunately, but it is the most accurate for the model I'm testing. -KG // note, where possible, arrays will all be js arrays, but can be operated on with math.js functions. // any math.js matrix objects will be labelled in comments // initialize counters let nsteps, nfailed, nfevals, npds, ndecomps, nsolves; nsteps = nfailed = nfevals = npds = ndecomps = nsolves = 0; // parse inputs const htspan = tspan.slice(0).reverse().reduce((acc, val, idx, arr) => idx < (arr.length - 1) ? val - arr[idx + 1] : acc, 0); const ntspan = tspan.length; const t0 = tspan[0]; var next = 1; const tfinal = tspan.slice(-1)[0]; const neq = y0.length; const eqIdx = math.range(0, neq); // ensure y0 is a column vector y0 = math.resize(y0, [neq, 1]); // check time direction, 1 ->, -1 <- const tdir = math.sign(tfinal - t0); let f0 = ode(t0, y0); // f0 is js array like y0 // options - thresholds const rtol = options.RelTol ? options.RelTol : 1e-3; const atol = options.AbsTol ? options.AbsTol : 1e-6; const threshold = atol / rtol; // options - hmax const safehmax = 16.0 * Number.EPSILON * Math.max(Math.abs(t0), Math.abs(tfinal)); const defaulthmax = Math.max(0.1 * Math.abs(tfinal - t0), safehmax); const himax = options.MaxStep ? options.MaxStep : defaulthmax; const hmax = Math.min(Math.abs(tfinal - t0), Math.abs(himax)); // increment nfun evals nfevals++; // refinement of output const refine = options.Refine ? Math.max(1, options.Refine) : 1; var outputAt, S; if (ntspan > 2) { outputAt = 'RequestedPoints'; } else if (refine <= 1) { outputAt = 'SolverSteps'; } else { outputAt = 'RefinedSteps'; S = Array.from({ length: refine - 1 }, (_, i) => (i + 1) / refine); } // mass matrix const Mt = math.identity(neq); // handle non-negative cases const idxNonNegative = options.NonNegative ? options.NonNegative : []; const nonNegative = idxNonNegative.length !== 0; let thresholdNonNegative; if (nonNegative) { [ode, thresholdNonNegative] = odenonnegative(ode, y0, threshold, idxNonNegative); // update f0 f0 = ode(t0, y0); nfevals++; } // prepare for loop let t = math.clone(t0); let y = math.clone(y0); let Mtnew = math.clone(Mt); // max Kth order: 5 is max let maxk = options.MaxOrder ? (options.MaxOrder <= 5 ? options.MaxOrder : 5) : 5; // Initialize method parameters. const G = [1, 3 / 2, 11 / 6, 25 / 12, 137 / 60]; // row const alpha = [-37 / 200, -1 / 9, -0.0823, -0.0415, 0]; // column const invGa = math.dotDivide(1, math.dotMultiply(G, math.subtract(1, alpha))).valueOf(); const erconst = math.add( math.dotMultiply(alpha, G), math.dotDivide( 1, math.resize(math.range(2, 6, true), [5, 1]) ) ).valueOf(); let difU = [ [-1, -2, -3, -4, -5], [0, 1, 3, 6, 10], [0, 0, -1, -4, -10], [0, 0, 0, 1, 5], [0, 0, 0, 0, -1] ]; const maxK = Array.from({ length: maxk }, (_, i) => i); const kJ = maxK.map(() => [...maxK]); const kI = maxK.map((v, i) => Array(neq).fill(v)); if (maxk < 5) { difU = difU.map(row => row.slice(0, maxk)).slice(0, maxk); } const maxit = 4; // get the jacobian matrix const jacThresh = math.add(math.zeros(math.size(y0)), atol); const [dfdy, nF] = odenumjac(ode, t, y, f0, jacThresh); // update fun evals counter nfevals += nF; // update npds npds++; // hmin is a small number such that t + hmin is clearly different from t in // the working precision, but with this definition, it is 0 if t = 0. let hmin = 16 * Number.EPSILON * Math.abs(t); //note: htry is empty & normcontrol false const wt = y.map(v => [Math.max(Math.abs(v), threshold)]); let rh = math.norm(f0.map((v, i) => [v[0] / wt[i][0]]), 'inf') * 1.25 / Math.sqrt(rtol); let absh = math.min(hmax, htspan); if (absh * rh > 1) { absh = 1 / rh; } absh = Math.max(absh, hmin); let h = tdir * absh; const tdel = (t + tdir * Math.min(Math.sqrt(Number.EPSILON) * Math.max(Math.abs(t), Math.abs(t + h)), absh)) - t; const f1 = ode(t + tdel, y); nfevals++; const dfdt = f1.map((v, i) => (v[0] - f0[i][0]) / tdel); const DfDt = math.multiply(dfdy, f0).map((v, i) => [v[0] + dfdt[i]]); // normcontrol => false, Mtype <= 0 rh = 1.25 * Math.sqrt(0.5 * math.norm(math.dotDivide(DfDt, wt), 'inf') / rtol); absh = Math.min(hmax, htspan); if (absh * rh > 1) { absh = 1 / rh; } absh = Math.max(absh, hmin); h = tdir * absh; // Initialize for main loop let k = 0; // first order bdf1 let K = 0; let klast = k; let abshlast = absh; let dif = math.zeros(neq, maxk + 2).subset(math.index(eqIdx, 0), math.multiply(h, f0)).valueOf(); let hinvGak = h * invGa[k]; let nconhk = 0; // step with current h,k counter let Miter = math.subtract(Mt, math.multiply(hinvGak, dfdy)).valueOf(); let LUp = math.lup(Miter); ndecomps++; let havrate = false; // Allocate output memory let chunk = Math.min(Math.max(100, 50 * refine), refine + Math.floor((2 ** 13) / neq)); let tout = ntspan > 2 ? math.zeros(1, ntspan).valueOf() : math.zeros(1, chunk).valueOf(); let yout = ntspan > 2 ? math.zeros(neq, ntspan).valueOf() : math.zeros(neq, chunk).valueOf(); let nout = 1; tout = math.subset(tout, math.index(0, nout - 1), t); yout = math.subset(yout, math.index(eqIdx, nout - 1), y); // Main Loop // init updated vars var psi, tnew, ynew, difkp1, invwt, minnrm, rate, oldnrm, err, nout_new, tout_new, yout_new; var done = false; var at_hmin = false; // let devcount = 0; // dev while (!done) { // ++devcount; // console.log(`::Main Loop:: (${++devcount})`); // dev hmin = 16 * Number.MIN_VALUE; absh = Math.min(hmax, Math.max(hmin, absh)); if (absh === hmin) { console.log('update: line 218'); return null; } else { at_hmin = false; } h = tdir * absh; // Stretch the step if within 10% of tfinal-t. if ((1.1 * absh) >= Math.abs(tfinal - t)) { h = tfinal - t; absh = Math.abs(h); done = true; } if ((absh !== abshlast) || (k !== klast)) { let difRU = math.multiply( cumprod( math.chain(kI) .subtract(math.multiply(math.add(kJ, 1), absh / abshlast)) .dotDivide(math.add(kI, 1)) .done() ), difU ); dif = math.subset( dif, math.index(eqIdx, K), // = math.multiply( math.subset(dif, math.index(math.range(0, neq), K)), math.subset(difRU, math.index(K, K)) ) ); hinvGak = h * invGa[k]; nconhk = 0; Miter = math.subtract(Mt, math.multiply(hinvGak, dfdy)); LUp = math.lup(Miter); ndecomps++; havrate = false; } // LOOP FOR ADVANCING ONE STEP. let nofailed = true // no failed attempts yet let wtcnt = 0; //dev while (true) { wtcnt++; // dev let gotynew = false; let ngoty = 0; // dev while (!gotynew) { ngoty++; // dev // Compute the constant terms in the equation for ynew. psi = math.resize( math.multiply( math.subset(dif, math.index(eqIdx, K)), math.multiply( math.subset(G, math.index(K)), math.subset(invGa, math.index(k)) ) ), [neq, 1] ); // predict new solution at t+h tnew = done ? tfinal : t + h; // purify h h = tnew - t; let pred = math.add(y, math.reshape(math.sum(math.subset(dif, math.index(eqIdx, K)), 1), [neq, 1])); ynew = math.clone(pred); // The difference, difkp1, between pred and the final accepted // ynew is equal to the backward difference of ynew of order // k+1. Initialize to zero for the iteration to compute ynew. difkp1 = math.zeros(neq, 1).valueOf(); invwt = math.dotDivide( 1, y.map((v, i) => [Math.max(Math.max(Math.abs(v[0]), Math.abs(ynew[i][0])), threshold)]) ); minnrm = 100 * Number.EPSILON * math.norm(math.dotMultiply(ynew, invwt), 'inf'); // Iterate with simplified Newton method for (var iter = 1; iter <= maxit; iter++) { let rhs = math.subtract( math.multiply(hinvGak, ode(tnew, ynew)), math.multiply(Mtnew, math.add(psi, difkp1)) ).valueOf(); // use math.lusolve to solve linear system let del = math.lusolve(LUp, rhs).valueOf(); let newnrm = math.norm(math.dotMultiply(del, invwt), 'inf'); difkp1 = math.add(difkp1, del); ynew = math.add(pred, difkp1).valueOf(); if (newnrm <= minnrm) { gotynew = true; break } else if (iter === 1) { if (havrate) { let errit = newnrm * rate / (1 - rate); if (errit <= 0.05 * rtol) { gotynew = true; break } } else { rate = 0; } } else if (newnrm > 0.9 * oldnrm) { break } else { rate = Math.max(0.9 * rate, newnrm / oldnrm); havrate = true; let errit = newnrm * rate / (1 - rate); if (errit <= 0.5 * rtol) { gotynew = true; break } else if (iter === maxit) { break } else if (0.5 * rtol < errit * rate ** (maxit - iter)) { break } } // update nrm oldnrm = newnrm; } // end of newton loop nfevals += iter; nsolves += iter; } // end while !gotynew //difkp1 is now the backward difference of ynew of order k+1. // MATLAB linear indexing could occur if k>maxk, though kmax let errColumn = Math.floor(k / maxk); err = math.multiply( math.norm(math.dotMultiply(difkp1, invwt), 'inf'), math.subset(erconst, math.index(k % maxk, errColumn)) ); if (nonNegative && (err <= rtol) && math.subset(ynew, math.index(idxNonNegative, 0)).some(v => v[0] < 0)) { console.log('Update line 339'); return null } if (err > rtol) { // % Failed step nfailed++; if (absh <= hmin) { return { t: math.subset(tout, math.index(0, math.range(0, nout))), sol: math.subset(yout, math.index(eqIdx, math.range(0, nout))), stats: { nsteps: nsteps, nfailed: nfailed, nfevals: nfevals, npds: npds, ndecomps: ndecomps, nsolves: nsolves } } } abshlast = absh; if (nofailed) { nofailed = false; let hopt = absh * Math.max(0.1, 0.833 * Math.pow((rtol / err), (1 / (k + 2)))); // 1/1.2 if (k > 0) { let errkm1 = math.chain(dif) .subset(math.index(math.range(0, math.size(dif)[0]), k)) .add(difkp1) .dotMultiply(invwt) .norm('inf') .multiply(erconst[k - 1]) .done()[0]; let hkm1 = absh * Math.max(0.1, 0.769 * Math.pow((rtol / errkm1), (1 / (k + 1)))); // 1/1.3 if (hkm1 > hopt) { hopt = Math.min(absh, hkm1); // don't allow step size increase k -= 1; K = math.range(0, k, true).valueOf(); } } absh = Math.max(hmin, hopt); } else { absh = Math.max(hmin, 0.5 * absh); } h = tdir * absh; if (absh < abshlast) { done = false; } let difRU = math.multiply( cumprod( math.chain(kI) .subtract(math.multiply(math.add(kJ, 1), absh / abshlast)) .dotDivide(math.add(kI, 1)) .done() ), difU ); dif = math.subset( dif, math.index(math.range(0, neq), K), math.multiply( math.subset(dif, math.index(math.range(0, neq), K)), math.subset(difRU, math.index(K, K)) ) ); hinvGak = h * invGa[k]; nconhk = 0; Miter = math.subtract(Mt, math.multiply(hinvGak, dfdy)); LUp = math.lup(Miter); ndecomps += 1; havrate = false; } else { // successful step break } } // end while true nsteps++; dif = math.subset(dif, math.index(eqIdx, k + 2), math.subtract(difkp1, math.subset(dif, math.index(eqIdx, k + 1)))); dif = math.subset(dif, math.index(eqIdx, k + 1), difkp1); for (let j = k; j >= 0; j--) { dif = math.subset(dif, math.index(eqIdx, j), math.add(math.column(dif, j), math.column(dif, j + 1))); } let NNreset_dif = false; if (nonNegative && math.subset(ynew, math.index(idxNonNegative, 0)).some(v => v[0] < 0)) { console.log("Update Line 421"); return null } // resolve this step for output switch (outputAt) { case 'SolverSteps': // computed points, no refinement nout_new = 1; tout_new = [tnew]; yout_new = [ynew]; break; case 'RefinedSteps': // computed points, with refinement let tref = S.map(v => t + (tnew - t) * v); nout_new = refine; tout_new = [...tref, tnew]; yout_new = [...ntrp15s(tref, tnew, ynew, h, dif, k, idxNonNegative), ynew]; break; case 'RequestedPoints': // output only at tspan points nout_new = 0; tout_new = []; yout_new = []; while (next < ntspan) { if ((tdir * (tnew - tspan[next])) < 0) { break; } nout_new += 1; tout_new.push(tspan[next]); if (tspan[next] === tnew) { yout_new.push(ynew); } else { yout_new.push(ntrp15s(tspan[next], tnew, ynew, h, dif, k, idxNonNegative)); } next += 1; } break; // switch break } if (nout_new > 0) { let oldnout = nout; nout = nout + nout_new; if (nout > tout.length) { tout = math.concat(tout, math.zeros(1, chunk).valueOf()); yout = math.concat(yout, math.zeros(neq, chunk).valueOf()); } let idx = math.range(oldnout, nout).valueOf(); tout = math.subset(tout, math.index(0, idx), idx.length > 1 ? tout_new : tout_new[0]); yout = math.subset(yout, math.index(eqIdx, idx), idx.length > 1 ? yout_new : yout_new[0]); } if (done) { break } // update step klast = k; abshlast = absh; nconhk = Math.min(nconhk + 1, maxk + 2); if (nconhk > k + 2) { var temp = 1.2 * (err / rtol) ** (1 / (k + 2)); var hopt = temp > 0.1 ? absh / temp : 10 * absh; var kopt = k; if (k > 0) { let errColumn = Math.floor((k - 1) / maxk); let errkm1 = math.multiply( math.norm(math.dotMultiply(math.subset(dif, math.index(eqIdx, k)), invwt), 'inf'), math.subset(erconst, math.index((k - 1) % maxk, errColumn)) ); temp = 1.3 * (errkm1 / rtol) ** (1 / (k + 2)); let hkm1 = temp > 0.1 ? absh / temp : 10 * absh; if (hkm1 > hopt) { hopt = hkm1; kopt = k - 1; } } // end k>0 if (k < maxk) { let errColumn = Math.floor((k + 1) / maxk); let errkp1 = math.multiply( math.norm(math.dotMultiply(math.subset(dif, math.index(eqIdx, k + 2)), invwt), 'inf'), math.subset(erconst, math.index((k + 1) % maxk, errColumn)) ); temp = 1.4 * (errkp1 / rtol) ** (1 / (k + 3)); let hkp1 = temp > 0.1 ? absh / temp : 10 * absh; if (hkp1 > hopt) { hopt = hkp1; kopt = k + 1; } } if (hopt > absh) { absh = hopt; if (k !== kopt) { k = kopt; K = math.range(0, k + 1).valueOf(); } } } // Advance integration 1 step t = tnew; y = ynew; if (NNreset_dif) { console.log("Update line 546"); return null } }// end main // setup output object: return { t: math.subset(tout, math.index(0, math.range(0, nout))), sol: math.subset(yout, math.index(eqIdx, math.range(0, nout))), stats: { nsteps: nsteps, nfailed: nfailed, nfevals: nfevals, npds: npds, ndecomps: ndecomps, nsolves: nsolves } } }34:Tab8,{ const time = solution.t; const sol = solution.sol; const sol2 = chunk.sol; // Parameters for the Model const colArea = typeParams.colArea; const caDark = typeParams.concCaDark; const cgDark = typeParams.concCgDark; const betaDark = typeParams.betaDark; const betaSub = typeParams.betaSub; const muRa = typeParams.muRa; const muTa = typeParams.muTa; const muPa = typeParams.muPa; const xi = typeParams.xi; // Cyclase parameters const nAlpha = typeParams.nAlpha; const rAlpha = typeParams.rAlpha; const KAlpha = typeParams.KAlpha / caDark; // Exchanger and Channel properties const fChCa = typeParams.fChCa; const nCh = typeParams.nCh; const KEx = typeParams.KEx / caDark; const KCh = typeParams.KCh / cgDark; // Rate for the Calcium dynamics const muCa = typeParams.muCa; const variables = ['Ra', 'Ta', 'Pa', 'y', 'z']; const ind_Ra = variables.indexOf('Ra'); const ind_Ta = variables.indexOf('Ta'); const ind_Pa = variables.indexOf('Pa'); const ind_y = variables.indexOf('y'); const ind_z = variables.indexOf('z'); // % cGMP and Calcium normalized by the dark values let cgmp = sol[ind_y].map(v => Math.exp(-v)); let ca = sol[ind_z].map(v => Math.exp(-v)); let cgmp2 = sol2[ind_y].map(v => Math.exp(-v)); let ca2 = sol2[ind_z].map(v => Math.exp(-v)); // % Currents normalized by the dark values // % intracellular recording (whole-cell patch) // % output.curr = 1-output.current; //let currCh = 2/(fChCa+2)*(1 - (1+KCh^nCh)./(cgmp.^nCh+KCh^nCh).*cgmp.^nCh); let currCh = cgmp.map(val => 2/(fChCa+2)*(1-(1+KCh**nCh)/(val**nCh + KCh**nCh)*val**nCh)); let currCh2 = cgmp2.map(val => 2/(fChCa+2)*(1-(1+KCh**nCh)/(val**nCh + KCh**nCh)*val**nCh)); //let currEx = fChCa/(fChCa+2)*(1 - (1+KEx)./(ca+KEx).*ca); let currEx = ca.map(val => fChCa / (fChCa + 2) * (1 - (1 + KEx) / (val + KEx) * val)); let currEx2 = ca2.map(val => fChCa / (fChCa + 2) * (1 - (1 + KEx) / (val + KEx) * val)); //let photocurrent = currCh + currEx; let photocurrent = currEx.map((v,i) => (v+currCh[i]) * maxResponse); let pc2 = currEx2.map((v,i) => (v+currCh2[i]) * maxResponse); // maxResponse is in pA let d = time[0].map((t,i) => {return {time:t,response:photocurrent[i]}}); let c = chunk.t[0].map((t,i) => {return {time:t,response:pc2[i]}}); return Plot.plot({ x: { label: "Time (s)" }, y: { label: "Photocurrent (pA)", grid: true }, marginTop: 30, marginBottom: 45, marks: [ Plot.ruleY([0]), Plot.ruleX([0]), Plot.lineY(d, {x: "time", y: "response"}), Plot.lineY(c, {x:"time",y:"response",stroke:"red"}) ] }); }35:Tb28,mdl = { const time = typeParams.time; const flashIntensities = typeParams.flashIntensities; const flashTimes = typeParams.flashTimes; const suppress_time = time .map(tVal => 1-(1-(Math.tanh((tVal-flashTimes[0])/0.0002) - Math.tanh((tVal-flashTimes[1])/0.0002))/2)) .map((v,i) => (time[i] < flashTimes[0]) ? 0 : v); const lightStimulus = Array(time.length) .fill(0) .map((val, i) => val + typeParams.flashIntensity * suppress_time[i]); // Parameters for the Model const colArea = typeParams.colArea; const caDark = typeParams.concCaDark; const cgDark = typeParams.concCgDark; const betaDark = typeParams.betaDark; const betaSub = typeParams.betaSub; const muRa = typeParams.muRa; const muTa = typeParams.muTa; const muPa = typeParams.muPa; const xi = typeParams.xi; // Cyclase parameters const nAlpha = typeParams.nAlpha; const rAlpha = typeParams.rAlpha; const KAlpha = typeParams.KAlpha / caDark; // Exchanger and Channel properties const fChCa = typeParams.fChCa; const nCh = typeParams.nCh; const KEx = typeParams.KEx / caDark; const KCh = typeParams.KCh / cgDark; // Rate for the Calcium dynamics const muCa = typeParams.muCa; const variables = ['Ra', 'Ta', 'Pa', 'y', 'z']; const varCount = variables.length; const ind_Ra = variables.indexOf('Ra'); const ind_Ta = variables.indexOf('Ta'); const ind_Pa = variables.indexOf('Pa'); const ind_y = variables.indexOf('y'); const ind_z = variables.indexOf('z'); // Initial conditions const initvalues = Array(varCount).fill(0).map(v => [v]); // Setup model function stimAmplitude(t) { return lightStimulus[time.reduce((sum, curr, i) => curr <= t ? i : sum, 0)]; } function odeFunc(t, u, ...odeArgs) { try { u = u.valueOf(); } finally { u = u.map(Number); } const ca = Math.exp(-u[ind_z]); const alpha = (1+Math.pow(KAlpha,nAlpha))/(rAlpha+Math.pow(KAlpha,nAlpha))*(rAlpha*Math.pow(ca,nAlpha)+Math.pow(KAlpha,nAlpha))/(Math.pow(ca,nAlpha)+Math.pow(KAlpha,nAlpha)); const dudt = new Array(u.length).fill(0); // Assuming '5' is the length of your 'u' array dudt[ind_Ra] = muRa * (colArea * xi * stimAmplitude(t) - u[ind_Ra]); dudt[ind_Ta] = muTa * (u[ind_Ra] - u[ind_Ta]); dudt[ind_Pa] = muPa * (u[ind_Ta] - u[ind_Pa]); dudt[ind_y] = u[ind_Pa] - betaDark * (Math.exp(u[ind_y]) * alpha - 1); dudt[ind_z] = muCa * ( ca * (1+KEx) / (ca+KEx) - (1+Math.pow(KCh,nCh)) / (Math.exp(-nCh * u[ind_y]) + Math.pow(KCh,nCh)) * Math.exp(-nCh * u[ind_y]) ) * Math.exp(u[ind_z]); return math.resize(dudt,[5,1]) } // ode options const options = odeoptions({'NonNegative':[0,1,2],'MaxStep':0.0002}); return {t:time,stim:lightStimulus,opt:options,odefun:odeFunc,init:initvalues} }36:[] 2d:{"title":"Photoreceptor Light Response Simulation","cells":[{"id":0,"value":"# Photoreceptor Light Response Simulation\nThis example simulates a light response in mouse rod and cone photoreceptors.","pinned":false,"hidden":false,"mode":"md"},{"id":2279,"value":"$2f","pinned":false,"hidden":false,"mode":"ojs"},{"id":1610,"value":"{\n let d = mdl.t.map((v,i)=>{return {time: v, stim: mdl.stim[i]};});\n return Plot.plot({\n marks: [\n Plot.ruleY([0]),\n Plot.lineY(d, {x: V=>V.time, y: V=>V.stim, height:200}),\n ]\n })\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":46,"value":"Rstar = 10;","pinned":false,"hidden":false,"mode":"ojs"},{"id":68,"value":"viewof cellType = Inputs.select([\"Rod\",\"Cone\"],{value:\"Rod\",label:\"Cell Type\"});","pinned":false,"hidden":false,"mode":"ojs"},{"id":1559,"value":"viewof maxResponse = {\n const inp = Inputs.range(\n [12,35],\n {\n value: 15,\n label:html`Maximal Response Amplitude (pA)`,\n step: 0.25,\n transform: x => x.toPrecision(2)\n }\n );\n return inp\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":11,"value":"## Imports","pinned":false,"hidden":false,"mode":"md"},{"id":180,"value":"math = await require(\"mathjs@11.8\");","pinned":false,"hidden":false,"mode":"ojs"},{"id":48,"value":"## Backend","pinned":false,"hidden":false,"mode":"md"},{"id":106,"value":"### Construct Stimulus Vector","pinned":false,"hidden":false,"mode":"md"},{"id":56,"value":"stimulusDuration = 0.01; // 10 msec","pinned":false,"hidden":false,"mode":"ojs"},{"id":108,"value":"responseDuration = 1.2; // sec","pinned":false,"hidden":false,"mode":"ojs"},{"id":111,"value":"sampleRate = 1000; // sampling rate","pinned":false,"hidden":false,"mode":"ojs"},{"id":114,"value":"stimulusOffset = 0.1;","pinned":false,"hidden":false,"mode":"ojs"},{"id":94,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":83,"value":"","pinned":false,"hidden":false,"mode":"html"},{"id":157,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":241,"value":"function odefinalize(statvect, nout, tout, yout) {\n let solver_output = {};\n solver_output.sol = tout; // Column vector\n solver_output.kvec = yout.valueOf();\n solver_output.dif3d = statvect; // Spread operator to convert array-like object to array\n return solver_output;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":255,"value":"function InfinityNorm(vector) {\n return Math.max(...vector.map(Math.abs));\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":247,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":261,"value":"function odenonnegative(ode, y0, threshold, idxNonNegative) {\n // Jacek Kierzenka\n // Copyright 1984-2020 The MathWorks, Inc.\n\n let thresholdNonNegative; \n if (threshold.length === 1) {\n thresholdNonNegative = idxNonNegative.map(() => threshold[0]);\n } else {\n thresholdNonNegative = idxNonNegative.map(index => threshold[index]);\n }\n \n function local_odeFcn_nonnegative(t, y, ...args) {\n let yp = ode(t, y, ...args);\n let ndx = idxNonNegative.filter(index => y[index] <= 0);\n ndx.forEach(index => {\n math.subset(yp,math.index(index,0),Math.max(math.subset(yp,math.index(index,0)), 0));\n });\n return yp;\n }\n return [local_odeFcn_nonnegative, thresholdNonNegative];\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1068,"value":"function cumprod(arr,dim=0) {\n // requires math.js\n let isMatrix = math.typeOf(arr) === 'DenseMatrix';\n arr = isMatrix ? arr.toArray() : arr;\n if (!dim) {\n arr = math.transpose(arr); // Transpose array to operate along columns\n }\n\n let result = arr.map((row) => {\n return row.reduce( (cump, val) => cump.concat((cump[cump.length-1] || 1) * val) , [] );\n });\n\n if (!dim) {\n result = math.transpose(result); // Transpose result back to original orientation\n }\n\n return isMatrix ? math.matrix(result) : result; // Convert back to matrix if original input was a matrix\n} ","pinned":false,"hidden":false,"mode":"ojs"},{"id":327,"value":"function odeoptions({\n AbsTol = null,\n BDF = null,\n Events = null,\n InitialStep = null,\n Jacobian = null,\n JConstant = null,\n JPattern = null,\n Mass = null,\n MassSingular = null,\n MaxOrder = null,\n MaxStep = null,\n NonNegative = null,\n NormControl = null,\n OutputFcn = null,\n OutputSel = null,\n Refine = null,\n RelTol = null,\n Stats = null,\n Vectorized = null,\n MStateDependence = null,\n MvPattern = null,\n InitialSlope = null\n} = {}) {\n return {\n AbsTol,\n BDF,\n Events,\n InitialStep,\n Jacobian,\n JConstant,\n JPattern,\n Mass,\n MassSingular,\n MaxOrder,\n MaxStep,\n NonNegative,\n NormControl,\n OutputFcn,\n OutputSel,\n Refine,\n RelTol,\n Stats,\n Vectorized,\n MStateDependence,\n MvPattern,\n InitialSlope\n };\n}\n","pinned":false,"hidden":false,"mode":"ojs"},{"id":1626,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2028,"value":"solution = await solve(mdl.odefun,mdl.t,mdl.init,mdl.opt);","pinned":true,"hidden":false,"mode":"ojs"},{"id":1827,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":2050,"value":"chunk = {\n const sol = solution.sol;\n const time = solution.t[0];\n const tspan = [0.15,0.4];\n const startIdx = time.reduce((sum, curr, i) => curr <= tspan[0] ? i : sum, 0);\n const endIdx = time.reduce((sum,curr,i) => curr <= tspan[1] ? i : sum, 0);\n const T1 = math.subset(time,math.index(math.range(startIdx,endIdx,true)));\n const y1 = math.subset(sol,math.index(math.range(0,5),startIdx));\n const chunk = await solve(mdl.odefun,T1,y1,mdl.opt);\n return chunk\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":288,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":2202,"value":"async function getChunk(soln,tspan){\n tspan = tspan.length === 2 ? tspan : [0,...tspan];\n const sol = soln.sol;\n const time = soln.t[0];\n const startIdx = time.reduce((sum, curr, i) => curr <= tspan[0] ? i : sum, 0);\n const endIdx = time.reduce((sum,curr,i) => curr <= tspan[1] ? i : sum, 0);\n const T1 = math.subset(time,math.index(math.range(startIdx,endIdx,true)));\n const y1 = math.subset(sol,math.index(math.range(0,5),startIdx));\n const chunk = await solve(mdl.odefun,T1,y1,mdl.opt);\n return chunk\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":14,"value":"## References","pinned":false,"hidden":false,"mode":"md"}],"files":"$W36","stdlib":"1","titling":"auto"}