Today I came across a good (or bad, depending on how you look at it) example of initial conditions / steady-states having a massive effect on model behaviour.

I think it is generally a good idea to get a model into a steady state* to use as its initial conditions (roughly speaking – as always – it never really reaches the true limit. To define steady state, we’ve used an L^2 norm of the difference in the state variable vector between subsequent paces, and said steady state is reached when this is < 1e-6, since you asked) .

Anyway, here’s a few reasons why:

• Physiological conditions are usually on steady-state timescales (years after the heart starts beating).
• Numerically speaking, you want a model that’s ‘stable’ and will keep going no matter how long the simulation is, so you want it to be capable of reaching a steady state.
• You have to start from somewhere, so steady state for some pacing rate has to be better than arbitrary numbers.
• [If you also simulate your experiment ‘to steady state’ then it looses its dependence on the initial conditions, so this would be a good thing if it’s a realistic portrayal of the experiment – if not then just simulate for the same amount of time as the experiment lasted.]
• If you aren’t in a steady state, the variables will drift off during your experiment, regardless of the experiment,  potentially influencing your results.

So I think it usually makes sense to get the model into a steady state before you do anything to avoid confusing ‘progress to steady state’ with what you are trying to simulate. Historically this wasn’t possible, when it took two weeks to compute an action potential; you had to put in your best guesses at initial conditions and just have a go. But nowadays (and certainly since the early 90s) that’s not really an excuse!

The model in question here is Priebe 1998 (CellML, Paper), one of the first human ventricular action potential models. I wanted to use it to demonstrate some of the capabilities of the Functional Curation Web Lab.

Now a lot of our virtual experiments use the model’s state variables taken when they’ve reached a steady state*, when pacing at a certain frequency, as the initial conditions for the reasons discussed above. We also tend to apply an intervention and then look at the steady state response. There were quite a few discrepancies with the results published in the original paper. For instance: the APD90 in our results was much longer, the S1-S2 curve correspondingly higher, and the result of IKr block completely different (all this initially got me worried about bugs in the CellML file, as we saw in the Decker model).

Figure 1 shows what it did under 1Hz pacing and IKr block on the Web Lab.

Figure 1: IKr block in the Priebe model when run for thousands of paces. Note that at 0% block (control) the APD90 is around 450ms. At high degrees of block the model fails to re-polarize. Have a look at these results on the web lab.

Figure 2 shows what it did when I just took the CellML file and did 10 paces (as they did in the paper – not steady state), and compared that with the results in the paper.

Figure 2: Left: my simulation with the same CellML file, just paced 10 times to control, then once per degree of drug block. Right: the result of the IKr block experiment in the original paper. Note the control APD90 is around 370ms. These plots look very similar. (click to enlarge)

Since the initial conditions don’t feature in the paper, it is impossible to tell for sure. But I think it’s safe to conclude that the paper contains results that are ‘on the way to steady state’. I’m not sure what drifts in the model to push the results elsewhere, it is usually an ionic concentration when it is a slow time scale thing.

There’s a slight subtlety here, in that you could say you were interested in the ‘immediate’ drug effect (Figure 1 shows the steady state drug effect, as well as using steady state initial conditions). In which case Figure 2 should look like this shouldn’t it? Well perhaps, but here the control case has an APD90 of ~370ms rather than the steady state ~440ms, so we can’t really say whether it would or not!

So my conclusions are:

• Think about what initial conditions to use:
  • if not steady state, have a good reason for this!
• Tell us what the initial conditions are (using steady state makes them easy to define!).
• Think about simulating ‘to steady state’:
  • if not, have a good reason for this!