Verification and testing of cellular ODE-based models

[michelebucelli]

Problem description

Both ionic models and active stress models are based on solving the same system of ODEs at every point in the domain. These systems are usually relatively complex, due to having many equations and involved expressions, so that making typos in their implementation is not uncommon.

As of now, the correctness of their implementation is only tested through regression tests, which verify that the solution of a 3D test case is consistent with a reference solution that was computed with the same code at the time of first implementation. This, however, has a few problems, in my opinion:

1. the comparison against the reference solution does not exclude that the reference solution is itself wrong (see e.g. Sign of stimulus current with Aliev-Panfilov model #535, where a typo in the implementation of the Aliev-Panfilov had gone unnoticed due to it being present when the reference solution was generated);
2. the regression tests usually involve less-than-trivial setups (e.g. the electrophysiology tests couple the ODE system with the monodomain PDE), so that if the test fails it is not immediately obvious whether the issue is in the ODE or in the PDE or elsewhere;
3. regression tests are usually 3D, which makes them relatively costly and slows down the CI/CD pipelines.

One alternative way of verifying the correctness would be to solve the ODE model as a standalone, and then compare its solution to that obtained by a trusted third-party implementation (cellular model developers often provide a reference implementation themselves). However, svMultiPhysics currently doesn't have a standalone ODE solver, as far as I can tell.

Proposed solution

I propose adding tests to the unitTests folder that implement a standalone 0D solver for both IonicModel and ActiveStress (or ActiveStressODE). These tests would take care of providing meaningful input to the ODE models surrogating coupling with 3D models (e.g. by solving an additional ODE to compute the membrane potential). The result would be compared against a reference solution provided e.g. as a CSV file (or other format).

The base classes IonicModel and ActiveStress should already support this thanks to the local timestepping functions (which might need to be moved to the public scope, or made available through friendship).

Additional context

1. @samibismar already did something of this kind as part of Add Regazzoni active stress model #596, although lacking the infrastructure that part was not committed anywhere. This proposed change would support more systematic testing.
2. I suspect that some of this overlaps with what @zasexton presented a while ago on auxiliary models (for the same reasons that IonicModel and ActiveStress have some overlap with it). I still think it might be useful to do this before auxiliary models come to the main branch, to provide testing infrastructure to cellular models that are currently under development, with the understanding that some of this might be discarded in favor of the more principled and general implementation.

[javijv4]

Hey @michelebucelli, I think this is a great idea! @kko27 mentioned @zasexton already has a similar issue in this; see here: #218

[michelebucelli]

@javijv4 thanks! I'd seen that issue before, and I think it would definitely be a useful addition! But its scope is much larger than what I'm proposing here, and I imagine that work on that will happen downstream of the quadrature refactoring (#593), assembly refactoring, and maybe others.

On the other hand, I think an automated testing framework for ODEs would be beneficial for currently ongoing development by me, @samibismar and @kko27 (and maybe others?), and could be implemented with relatively minimal additions, even if not in the most general way.

[javijv4]

@michelebucelli, I agree! I think it is a good idea to set this up now, since a lot of people are working on new ODE models.

[samibismar]

I have added isolated GTest here coverage for RegazzoniActiveStress, including initialization, comparison of the full 20-state trajectory and active tension against the patched Regazzoni reference implementation, and RU probability conservation/finiteness checks.
The Regazzoni-specific tests pass, and the full unit-test suite currently passes 236/236.
I’ll use this as the starting point for extending isolated verification to the other active-stress and ionic models.

[michelebucelli]

Thank you @samibismar! I had a look at the code, and here's a few suggestions.

1. Right now, you have a custom "harness" for every one of the two models, to expose some of their private methods so that the standalone ODE solver can call them. I think this scales sub-optimally as new models are added. What I suggest is to move the local time stepping methods from ActiveStress into the public scope of that class (I don't see any problem with that, since they do not change the class' state as far as I recall), and then write the ODE solver so that it uses polymorphism to be independent of the concrete model being solved.
2. For NashPanfilov, you use the harness to configure the model parameters as well. I think you can do that by changing the values of NashPanfilov::Parameters instead, before calling read_model_specific_parameters.
3. The reference solutions are now hardcoded into the test file. I suggest to read them from an external CSV file (or equivalent format), so that data and logic are cleanly separated and updating the reference solution will not require recompilation.

[samibismar]

@michelebucelli I addressed the suggestions from your review:

• removed the per-model active-stress test harnesses and now exercise the node-local methods polymorphically through ActiveStress&
• added a public parameter-configuration wrapper while keeping the model-specific implementation protected
• Nash-Panfilov parameters are now configured through NashPanfilov::Parameters
• moved the Nash-Panfilov and Regazzoni trajectory reference data into external CSV files
• added shared active-stress trajectory test infrastructure

All unit tests pass after the refactor (240/240).

I’m moving on to the ionic-model tests next. I wanted to check the preferred testing interface there. IonicModel already exposes init() and integ() publicly, but getf() remains protected. For direct RHS verification, would you prefer exposing the local/RHS interface similarly toActiveStress, or should the ionic tests exercise only the existing public integration interface?