2.3.1. Integrating systems of differential equations

In this tutorial we look into the integration of systems of differential equations, using example models from the domain of computational physiology. We also examine the the main control parameters that can have a significant impact on the performance and accuracy of a given method.

2.3.1.1. Numerical integrators

2.3.1.1.1. Euler method

The Euler method to integrate a system of ordinary differential equations is probably the most well known method, particularly popular given its simplicity. Due to its simplicity it is, however, usually not the most appropriate method to use when performing simulation experiments with complex models of biological systems. It is, however, a very useful example to use to demonstrate key concepts that apply to most numerical integration methods.

Given the initial value problem

(1)\[y'(t) = f(t, y(t)), \quad\quad y(t_0) = y_0,\]

we now want to approximate \(y(t)\) over some interval \(t_0 \rightarrow t_{final}\). The Euler method approximates \(y(t)\) by dividing the interval into a series of steps, of size \(h\), and stepping through the interval. One step of the Euler method from \(t_n\) to \(t_{n+1} = t_n + h\) is given by

(2)\[y_{n+1} = y_n + h f(t_n, y_n).\]

The value of \(y_n\) is an approximation of the solution of the ODE (1) at time \(t_n\): \(y_n \approx y(t_n)\). The Euler method proceeds through the interval in constant steps of size \(h\) until \(t_n = t_{final}\) and the integration is complete.

2.3.1.1.1.1. Task 1 - the effect of step size

As you’d imagine, the size of \(h\) in the Euler method is crucial to the successful application of the method. For the successful (although perhaps inaccurate) integration of a typical physiological model (see the Physiome Repository for a collection of examples), \(h\) can be so small that the computational cost of performing the simulation is very high. In some cases, mathematical models may be unsuitable for integration by Euler method, regardless of how small the step size is reduced to.

In the follow task, we investigate the effect of altering the size of \(h\) on two separate simulation experiments. The first looks at a simple mathematical model where the experiment could be replicated with pen and paper, while the second looks into the application of Euler’s method to a biophysical model of cellular electrophysiology.

In this task we use the trivial model:

(3)\[\begin{split}x(t) &= \mathit{sin}(t),\\ y'(t) &= \mathit{cos}(t) \quad\mathrm{with}\quad y(0) = 0.\end{split}\]

As you can see from (3), if correctly integrated \(x(t)\) and \(y(t)\) should be identical. We now use this model to demonstrate the effect of step size (\(h\)) on simulating this model using Euler integration over the interval \(0 \rightarrow 2 \pi\).

1. Run MAP Client, choose File ‣ Open and select HOME/projects/mapclient-workflows/DTP-Simulation-Task1.
2. This simple workflow should look similar to Fig. 2.16. The workflow is pre-configured so there is no configuration required.

Fig. 2.16 The first Euler example as loaded.

3. Click the Execute button and you should get a widget displayed as per Fig. 2.17.

Fig. 2.17 The cool Euler integrator interface. In this simple interface, you will see the standard sine function, \(sin(t)\), plotted in the right hand panel. The toolbar under the plot is self-explanatory, but provides access to some nifty features. At the top of the left hand panel you will see the control to set the Euler step size for this model, \(h\) and also the number of points to be obtained. The Simulate button will execute the Euler integration of the model (3) and plot the result on the plot to the right. This can be repeated with various values of \(h\). The Clear graph button will, surprisingly, clear the current simulation results from the plot panel. The Done button will drop you back to the work-flow diagram, where you can get back to the plot by executing the work-flow once more.

4. As described in Fig. 2.17, multiple simulations can be performed with varying values for the step size, \(h\). Shown in Fig. 2.18 you can see that as \(h\) reduces in size, the approximation of the model (3) by integration using the Euler method gets more accurate.

Fig. 2.18 Simulation results demonstrating the effect of step size, \(h\), on the accuracy of Euler’s method in approximating the solution of (3).

5. Now have a play with combining different values for the step size and the number of points to be obtained. See if you can answer the following.
   1. How small should \(h\) be to accurately simulate a sine wave?
   2. What do you think would happen beyond a single cycle?
   3. Given \(h = 1\), do you obtain a more accurate solution with a large number of points or a small number of points?

2.3.1.1.2. CVODE

From the Sundials suite of tools, CVODE is a solver for stiff and nonstiff ordinary differential equation (ODE) systems (initial value problem) given in explicit form in (1) above. CVODE is widely regarded as one of the gold standard implementations of a robust and flexible numerical integrator. One of the advantages of CVODE over Euler’s method is that it makes use of adaptive stepping over the interval of integration - rather than taking fixed sized steps through time, for example, CVODE will determine how quickly things are changing and adjust the size of the step accordingly.

2.3.1.1.2.1. Task 2 - fixed vs adaptive stepping

In this task we examine the limitations and the computational costs associated with a fixed step method (Euler) compared to an adaptive step method (CVODE). Here we continue with our sine integration demonstration model to help highlight the differences.

1. Run MAP Client, choose File ‣ Open and select HOME/projects/mapclient-workflows/DTP-Simulation-Task2.
2. This simple workflow should look similar to that used in task 1 above (Fig. 2.16). The workflow is pre-configured so there is no configuration required.
3. Click the Execute button and you should get a widget displayed as per Fig. 2.19.

Fig. 2.19 The user interface in this task is similar to that described in Fig. 2.17, and the common elements behave the same. In addition there is the ability to choose either the Euler or CVODE numerical integration methods. As the CVODE method is an adaptive stepping method, the value of \(h\) is used to limit the maximum step size that the algorithm will use, with \(h=0\) indicating the maximum step size is unlimited.

4. You can now easily see the difference between the two integration methods by directly comparing them, as shown in Fig. 2.20.

Fig. 2.20 Simulation results showing the comparison between the Euler and CVODE integrators.

5. Now have a play with step sizes, number of points, and integration methods to explore the features of these two integration methods and see if you can address these questions.
   1. What is the largest maximum step size you can use with CVODE to accurately simulate a sine wave with number of points being set to 2?
   2. How small does \(h\) need to be to get the same solution with Euler?
   3. Are either of those a useful solution?
   4. What is the minimum number of points required to capture an accurate sine wave?
   5. Can you determine a configuration for Euler and CVODE which demonstrates a cheaper, more accurate, simulation using CVODE with this model?

2.3.1.1.2.2. Task 3 - error control

In addition to providing adaptive stepping, CVODE is also a very configurable solver. Beyond the maximum step size explored above, a further control parameter of that is often of interest are the tolerances used to control the accumulation of error in the numerical approximation of the mathematical model. This tolerance specifies how accurate we require the solution of the integration to be, and the value used can be very specific to the mathematical model being simulated. In task 2 above, we used a tolerance of 1.0e-7. Depending on the behaviour of your mathematical model, you may need to tighten (reduce) or loosen (increase) the tolerance values, depending on the specific application the model is being used for. Here we explore the effect of the tolerance value on the ICC model introduced above.

We use the recent biophysically based mathematical model of unitary potential activity in interstitial cells of Cajal. The interstitial cells of Cajal (ICC) are the pacemaker cells of the gastrointestinal tract and provide the electrical stimulus required to activate the contraction of smooth muscle cells nescessary for the correct behaviour of the GI tract. This particular model was developed by scientists at the Auckland Bioengineering Institute to investigate a specific hypothesis regarding the biophysical mechanism underlying the pacemaker function of ICCs.

1. Run MAP Client, choose File ‣ Open and select HOME/projects/mapclient-workflows/DTP-Simulation-Task3.
2. This simple workflow should look similar to that used in task 1 above (Fig. 2.16). The workflow is pre-configured so there is no configuration required.
3. Click the Execute button and you should get a widget displayed as per Fig. 2.21.

Fig. 2.21 The user interface in this task is similar to that described in Fig. 2.17, and the common elements behave the same. We now are only using the CVODE integration method so \(h\) is the maximum step size with \(h=0\) indicating an unlimited step size. The tolerance value for the simulation can also be edited in this interface.

4. You can now investigate the effect of changing the tolerance value and maximum step size on the simulation result. Not all combinations will successfully complete. Example results are shown in Fig. 2.22.

Fig. 2.22 Simulation results for a selection of simulations of the ICC model using various configurations of the CVODE integratior.

5. After exploring the effects of the integrator parameters and the simulated model behaviour, see if you can answer the following questions.
   1. With \(h=0.0\), how loose can the tolerance be and still get an accurate solution?
   2. How tight can you make the tolerance before the computational cost outweighs any improvement in solution accuracy?
   3. Is there any value of \(h\) that will give an accurate solution for a tolerance of 0.01?