nmodl icon indicating copy to clipboard operation
nmodl copied to clipboard

Side effects in KINETIC blocks are delayed.

Open 1uc opened this issue 1 year ago • 0 comments

Let's start by reviewing Newtons method for F(X) = 0 with Jacobian J(X):

while(true) {
  F, J = compute_f_and_jacobian(X);
  if(is_converged(X, F)) { break; }

  // Rest of Newton to compute the
  // next `X`.
}

the important point: We run compute_f_and_jacobian once more after having found the converged values of the state variables.

Consider the following KINETIC block:

KINETIC state {
  ~ pump + pumpca <-> cao (1, 1)
  ~ ca <-> pumpca (1, 1)
  ica = f_flux - b_flux
  cai = ca
}

(loosely related to cabpump.mod.)

AFAICT, in NOCMODL the "side effects" are part of the compute_f_and_jacobian. Hence, they get run once more after the states have been updated.

In NMODL we inline the meaning of {f,b}_flux and then move all "side effects" upwards:

KINETIC state {
  i = /* expression involving ca, pumpca & rates */
  cai = ca

  ~ pump + pumpca <-> cao (1, 1)
  ~ ca <-> pumpca (1, 1)
}

Then we create a non-linear system from the contiguous block of reaction equations. We then put all the "side effects" in a method called initialize. The function compute_f_and_jacobian doesn't contain the side effects. The only assignments are:

F[...] = ...;
J[...] = ...;

We then do:

initialize();
newton(compute_f_and_jacobian, ...);
finalize(); // (always) empty

Which means that the side effects are computed only once, before the first iteration of Newton. Therefore, they lag behind by one timestep.

This happens for statements like:

cai = ca

but also:

ica = f_flux - b_flux

1uc avatar Aug 09 '24 09:08 1uc