modcc: ion concentration writers

[noraabiakar]

Note: the internal concentration of "ca" is used here, but the issue applies for all ions and for both internal and external concentrations

How do nrnivmodl and modcc interpret this nmodl code?

NEURON {
    SUFFIX ca_write
    USEION ca READ ica WRITE cai VALENCE 2
}
ASSIGNED { ica }
STATE    { concentration cai }
INITIAL  { concentration = cai }

nrnivmodl

• assigns concentration to the internal concentration of the "ca" ion.
• Produces similar code if cai is a PARAMETER or ASSIGNED instead of STATE, and if cai is READ or READ/WRITE instead.

modcc

• Requires a few changes: a BREAKPOINT block, and for ica to not be ASSIGNED so :

NEURON {
    SUFFIX ca_write
    USEION ca READ ica WRITE cai VALENCE 2
}
STATE    { concentration cai }
INITIAL  { concentration = cai }
BREAKPOINT {}

• Generates code that assigns concentration to a mechanism internal variable named cai, which is not going to contain the value of the internal concentration of "ca". The same code is generated if cai has READ instead of WRITE permission.

void init(mechanism_cpu_ca_write_pp_* pp) {
    int n_ = pp->width_;
    for (int i_ = 0; i_ < n_; ++i_) {
        pp->concentration[i_] = pp->cai[i_];
    }
}

• If cai is not a STATE variable (or anything else) and it has WRITE permission: generates code that expect the internal concentration to be written:

  void init(mechanism_cpu_ca_write_pp_* pp) {
  int n_ = pp->width_;
  for (int i_ = 0; i_ < n_; ++i_) {
      auto ion_ca_index_i_ = pp->ion_ca_index_[i_];
      ::arb::fvm_value_type cai = 0;
      pp->concentration[i_] = cai;
      pp->ion_ca_.internal_concentration[ion_ca_index_i_] = fma(pp->weight_[i_], cai, pp->ion_ca_.internal_concentration[ion_ca_index_i_]);
  }
  }

• If cai is not a STATE variable (or anything else) and it has READ permission: generates code that assigns concentration to the internal concentration of the "ca" ion.

  void init(mechanism_cpu_ca_write_pp_* pp) {
  int n_ = pp->width_;
  for (int i_ = 0; i_ < n_; ++i_) {
      auto ion_ca_index_i_ = pp->ion_ca_index_[i_];
      ::arb::fvm_value_type cai = pp->ion_ca_.internal_concentration[ion_ca_index_i_];
      pp->concentration[i_] = cai;
  }
  }

• cai can't have both READ and WRITE permissions in the nmodl code.

---

Summary of modcc's behavior:

• When cai has READ permission and cai is not a STATE variable, we observe the correct behavior: variables assigned to cai, get assigned the internal concentration of "ca".
• The presence of cai in STATE complicates things: cai no longer points to the internal concentration of "ca", but a separate uninitialized (or zero'd) variable. READ permission no longer works, and WRITE permission doesn't give read access either. (Question: what does cai being a STATE variable mean here? It alters the meaning of cai which now indicates it is both the internal concentration of "ca" and not at the same time.)
• When cai has WRITE permission and cai is a STATE variable, modcc expects cai to be initialized in the mechanism, the user can't expect it to have the same value as the internal concentration of "ca" if it is not initialized.
• When cai has WRITE permission and cai is not a STATE variable, modcc expects the internal concentration of "ca" to be written anytime it is used, which means that if it is read in some places but not written, the internal concentration will be written to zero.

Expected behavior of modcc:

• Our specifications of nmodl indicate that if a variable has WRITE permission, it has READ permission as well. It's reasonable for users to expect that variables assigned to cai get assigned the internal concentration of "ca".
• When cai has WRITE permission, it's unreasonable for the internal concentration of "ca" to be modified if cai is not.
• The classification of cai as a STATE variable should be better explained, and it should result in sane code for both READ and WRITE permissions, or raise an exception.

[noraabiakar]

modcc's current behavior has been designed to support common nmodl patterns:

• cai is usually either read or written, but not both.
• When cai is read, it is not usually a STATE variable.
• When cai is written, it is initialized within the mechanism, and doesn't rely on the initial internal concentration stored in the simulator.
• When cai is written, it is usually integrated, not directly assigned.

If these are the only patterns we want to support, we need to properly document this and ideally also raise exceptions if users attempt to do something different.

[thorstenhater]

Points to make clear in the documentation (at least for me, these were confusing. If I got some of them still wrong, all the more reason to document them) (also sorry for repeating some from above)

• There are two Xi in play, having the same name. One is the actual, cell state, one a shadow copy in the mechanism.
• If Xi is WRITE within the mechanism, this is added to the actual Xi
• If Xi is READ however, we will give you the actual value
• Multiple mechanisms might execute in arbitrary order and/or concurrently.
• Therefore, the pattern Xi += v (and variants thereof) is ill-defined and we rather have you not do it
• WRITE is not read-write. Reading Xi will give you the shadow value. Which might not be initialised prior to the first write.

Suggestions/Thoughts:

• Additive contributions should be zero-init'ed
• Suspiciuous behaviour like listed above should at least trigger a warning, possibly an error.
• Having Xi in STATE (also PARAMETER, ASSIGNED) seems dubious from a model standpoint, or am I missing something?
• Naive question: Do we have examples were nrnivmodl's behaviour is used to some (sensible) effect?
• Very naive question: I would assume people modelling exponential behaviour will instinctively write Xi = c Xi. However, that is better(?) done using Xi' = c, but that would require Xi being STATE. What is our answer to a request for modelling this?
• Personally, I would prefer to enable the pattern Xi = c Xi by differentiating between global and shadow copies. Something like this maybe

  NEURON {
  SUFFIX ca_write
  USEION ca READ cai AS cai_global WRITE cai AS cai_contribution VALENCE 2
  }

  Obvious drawback not mentioned.

[halfflat]

• If Xi is WRITE within the mechanism, this is added to the actual Xi
• If Xi is READ however, we will give you the actual value
• Multiple mechanisms might execute in arbitrary order and/or concurrently.

It's only added (with a weight factor) in order to handle partial CV coverage by a concentration-writing mechanism. We explicitly test to ensure that no two mechanisms writing to the same ionic concentration overlap (see fvm_layout.cpp lines 920-937).

This means that once we make it past the initialization phase, there should be no order-dependence issues in ion concentration writes.

Note that for each concentration, we compute in fvm_build_mechanism_data two values: an 'initial concentration', and a 'reset concentration'. In the integration loop, when it comes time to update the concentrations, we set it first to the 'initial' value, and then add the weighted contributions from the mechanisms. This initial value is the contribution from the concentration values painted on the cell (or from the defaults) in a CV, wherever there isn't some mechanism writing it. The reset value on the other hand is what the concentration would be without any mechanisms writing it all. This is used to set the CV concentration values before mechanism initialization (and after a call to reset()), so that mechanisms that try to read the concentration in their INITIAL block see something sensible.

• Additive contributions should be zero-init'ed

This is what we're doing with the process outlined above.

[halfflat]

Just to clarify the semantics as implemented: for Arbor, a concentration-writing mechanism doesn't 'see' the effects of any other concentration-writing process, or averaging, etc. It is responsible purely for governing the evolution of the concentration as a function of ionic currents and any internal state it might hold.

Consequently, the only way the mechanism should be able to read that ionic concentration is at mechanism initialization time, and the value it should see at initialization is the 'reset concentration' described above.

[thorstenhater]

The request for zeroing these contributions was motivated by the idea that a user -- knowing it is an additive component -- might not explicitly initialise Xi and read from it later on in the mechanism. C++ people might be aware that anything can happen here, but that might not be true of scientists writing models.

[halfflat]

Inasmuch as Xi is a state variable that can be read in some contexts, it should be explicitly initialized to something useful (like the external Xi init value). Without checking the code — yes, I know this is a sin — I thought all our mechanism state variables were initialized to NaN before dealing with the INITIAL block, so we shouldn't be seeing uninitialized values in the state Xi, just NaN values.

The issue is, in part, NMODL. When we see Xi in a mechanism, depending on context it can mean the external ion concentration, a random state variable or parameter unrelated to an ion, or a state variable that represents an ion concentration, from which the external ion concentration is taken. If we want concentration writers to be able to read the default ion concentration value too, then Xi can mean both the state and the external in different spots within the same mechanism. It's just ugly.

But again, we should determine what semantics we want and proceed accordingly.

[thorstenhater]

As we have found out (and documented) over implementing the #1729 feature, this is part of the cable model, if a different model is required, one can use the diffusive concentration Xd, without actual diffusion as needed.