Sign of stimulus current with Aliev-Panfilov model

[michelebucelli]

Description

As part of working on #530, I tried running the test tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_pfib_AP.xml over multiple timesteps. The results seem to be incorrect.

In the initial timesteps, the excited region's potential becomes lower than the resting value (instead of higher):

As the simulation proceeds, the (negative) excitation does not propagate, and at some point it blows up, likely due to the negative values falling in some unfeasible range for the model (next images are the same frame but with different color scales):

The simulation crashes in time step 42, throwing the error std::runtime_error: [cep_integ_l] A NaN has been computed during time integration of electrophysiology variables.

Reproduction

On the main branch, run the test case tests/cases/cep/niederer_benchmark_ECGs_quadrature/solver_GMRES_FE_pfib_AP.xml setting Number_of_time_steps to 1000 (or any value above 40).

Expected behavior

The excited region should have higher potential then the resting value, as in this picture (same timestep as first picture above):

The stimulus should propagate over time to the whole domain, and should not blow up.

The above was achieved by changing the sign of the applied current stimulus from negative to positive. Accordingly, I suggest that that sign be changed in the parameter file for this test and possibly other tests involving the AP model.

Additional context

I don't know if the same issue happens with models other than AP.

I noticed that the sign of the applied current is interpreted differently for different models (at least for TTP and AP, not sure about the others). I think it might be beneficial to make this setting consistent, so that the applied current parameter section can become independent of the ionic model being used. If people agree with this, I can take care of it as part of the refactoring proposed in #530.

Code of Conduct

• I agree to follow this project's Code of Conduct and Contributing Guidelines

[michelebucelli]

Partially related to this issue, I saw that the value of CepModAp::alpha, which corresponds to the parameter $a$ in the original Aliev-Panfilov paper, is set to $10^{-2}$, instead of $0.15$ as in the paper. Do you know if there's a reason for this change, and/or a different source paper for this change?

[alisonmarsden]

@vvedula22 ran this benchmark case previously using the Fortran solver and he might have input on this issue.

[vvedula22]

Both Aliev-Panfilov and Fitzhugh-Nagumo models were adapted from Ellen's work. They have an entire phase plot for various parameters, and I believe I have tested multiple parameters and eventually retained $\alpha = 0.01$. Here is a case that uses the same parameters as given in the original work by Goktepe with $\alpha = 0.05$.

Regarding Niederer's benchmark, you're right. The stimulus amplitude is treated differently, partly due to the way the equations were written in the original papers and when they were implemented and tested. For e.g., in Goktepe et al., the external stimulus is added directly as a source term ($f_{ext}$). However, in the works by ten-Tusscher et al. (2004, 2006), the stimulus current appears on the right-hand side with a negative sign (i.e., $$\frac{dV}{dt} = -\frac{1}{C_m} \left(I_{ion} + I_{stim}\right)$$). I have only tested Niederer's benchmark with the ten-Tusscher-Panfilov model, not others.

I believe apart from Aliev-Panfilov and Fitzhugh-Nagumo models, all other models (Bueno-Orovio, Nygren, Stewart, TTP) require a negative stimulus current.

[michelebucelli]

Thank you for the clarifications @vvedula!

What do you think about using the same sign for the applied current across all ionic models? It would amount to switching the sign in the formulation of the Aliev-Panfilov and FitzHugh-Nagumo models, so that all should work correctly with the same applied current value. This should come naturally as part of #530.

I would be inclined to use a positive sign (so that the user does not need to worry about it) and put the minus somewhere in the source code, but I'm open to leaving it negative for all models if that is more backwards-compatible.

[alisonmarsden]

I think they should all be consistent. Let's pick whichever convention is most commonly used in the literature for the models we plan to use the most (e.g. TTP). @MatteoSalvador might have comments here as well!

[vvedula22]

Generally, we implemented the Aliev-Panfilov and Fitzhugh-Nagumo models at the very beginning as we were learning the basics of electrophysiology, and to test the operator-splitting approach in our code. We never really used it later, but focused on the other models (TTP, BO, etc.).

I'd suggest you generalize the stimulus to this equation:

$$\frac{dV}{dt} = -\frac{1}{C_m} \left( I_{ion} + I_{stim} \right)$$

as this is the standard convention used in the literature. At the same time, in the documentation, it would be helpful to note that the stimulus convention for these two excitation models (AP, FN) is reversed from that published in the literature. You may also display the equation with a negative forcing term for these two models.

[alisonmarsden]

Sounds good to me. Tagging @kko27 to make sure he's aware!