Iterative Implicit Methods for Solving Hodgkin-Huxley Type Systems

TL;DR

This paper investigates implicit numerical methods for solving Hodgkin-Huxley type models, demonstrating their ability to preserve system dynamics and improve accuracy through adaptive time stepping, especially in chaotic regimes.

Contribution

It introduces implicit methods tailored for Hodgkin-Huxley models, analyzing their bifurcation behavior and effectiveness in capturing complex dynamics.

Findings

Implicit methods preserve limit cycles in the model.

Adaptive time stepping enhances accuracy in chaotic regions.

The approach aids in understanding bifurcation patterns.

Abstract

We are motivated to approximate solutions of a Hodgkin-Huxley type model with implicit methods. As a representative we chose a psychiatric disease model containing stable as well as chaotic cycling behaviour. We analyze the bifurcation pattern and show that some implicit methods help to preserve the limit cycles of such systems. Further, we applied adaptive time stepping for the solvers to boost the accuracy, allowing us a preliminary zoom into the chaotic area of the system.

Tables

Table 1. Table 1: Runtimes of the different algorithms for the interval [ 0 , 10000 ] 0 10000 [0,10000] . The computations were carried out on an Intel Core i5-7200U.

	S=100			S=258.15
Scheme	min time	mean time	max time	min time	mean time	max time
ISIE	296.279 ms	299.401 ms	307.805 ms	263.207 ms	266.759 ms	273.071 ms
ICN	467.571 ms	473.274 ms	479.486 ms	405.481 ms	415.413 ms	453.097 ms
ISV	428.761 ms	435.463 ms	450.706 ms	393.316 ms	396.468 ms	401.683 ms
MMRK4	434.206 ms	451.805 ms	472.120 ms	403.617 ms	414.823 ms	429.899 ms
IRK4	1.149 s	1.153 s	1.160 s	972.086 ms	1.065 s	1.135 s
PIDICN	365.669 ms	371.045 ms	382.652 ms	128.728 ms	144.721 ms	203.532 ms
AICN	929.919 ms	937.354 ms	954.432 ms	689.332 ms	743.694 ms	823.602 ms
ARK4	61.272 ms s	64.606 ms	70.739 ms	51.398 ms	60.307 ms	89.615 ms
AIRK4	392.592 ms	398.475 ms	413.562 ms	373.588 ms	397.572 ms	434.646 ms

Table 2. Table 2: A tabular view of the previously defined global norms for all solvers. The left side of the table contains the stable cycling case S = 100 𝑆 100 S=100 while the right side contains the chaotic cycling case S = 258.15 𝑆 258.15 S=258.15 .

	S=100		S=258.15
Scheme	${\parallel 𝐮\parallel }_{L_2[0,9999]}$	${\parallel \frac{d𝐮}{dt}\parallel }_{L_2[0,9999]}$	${\parallel 𝐮\parallel }_{L_2[0,9999]}$	${\parallel \frac{d𝐮}{dt}\parallel }_{L_2[0,9999]}$
ISIE	1779.1042482178611	763.6817995346803	2074.99633068202	551.5939982581531
ICN	1779.6596093842356	765.6633264308293	2080.081103110141	561.9538618698593
ISV	1779.8227558934245	764.772539199622	2072.520439246565	548.0642069057208
MMRK4	1779.4729101969149	764.9818721864283	2066.8111537452933	538.8410380757709
IRK4	1779.659878381908	765.6883397814581	2094.416156306047	584.595355425552
PIDICN	1778.7230770630867	765.3645009066148	2086.96979063115	620.090939713537
AICN	1779.6596093842356	765.6633264308293	2080.081103110141	561.9538618698593
ARK4	1772.6310255573324	764.9749934137184	2078.7562825863693	565.8865895392348
AIRK4	1773.8147243055737	765.2310905681921	2068.84774687211	546.2333427143521

Table 3. Table 3: A tabular view of the statistical features for the adaptive solvers. The left side of the table contains the stable cycling case S = 100 𝑆 100 S=100 while the right side contains the chaotic cycling case S = 258.15 𝑆 258.15 S=258.15 .

	S=100		S=258.15
Scheme	$𝔼\left[{\parallel 𝐮\parallel }_{L_2[0,9999]}\right]]$	$𝕍\left[{\parallel 𝐮\parallel }_{L_2[0,9999]}\right]$	$𝔼\left[{\parallel 𝐮\parallel }_{L_2[0,9999]}\right]$	$𝕍\left[{\parallel 𝐮\parallel }_{L_2[0,9999]}\right]$
PIDICN	1.53167	3.94941	1.56736	8.5948
AICN	1.10589	1.3043	0.934257	0.936426
ARK4	3.64771	14.896	2.74606	11.2035
AIRK4	3.33507	12.3004	2.60268	9.33504

Equations

$$C\frac{\partial V}{\partial t}$$

$$\frac{dn}{dt}$$

$$\frac{dm}{dt}$$

$$\frac{dh}{dt}$$

$$n_{\infty }(V)$$

$$m_{\infty }(V)$$

$$h_{\infty }(V)$$

$${\alpha }_n(V)$$

$${\alpha }_m(V)$$

$${\alpha }_h(V)$$

$$C\frac{dV}{dt}$$

$$\frac{dn}{dt}$$

$$\frac{dm}{dt}$$

$$\frac{dh}{dt}$$

$$\dot{V}$$

$$\dot{w}$$

$${\tau }_o(o)\frac{do}{dt}=$$

$${\tau }_i(o)\frac{d{a}_i}{dt}=$$

$${\tau }_o(o)\frac{do}{dt}=$$

$${\tau }_i(o)\frac{d{a}_i}{dt}=$$

$$\begin{array}{rl}{\tau }_x\frac{dx}{dt}& =-x-\sum _{i\in \{he,li,le\}}{a}_i{w}_i(x-x_i)-a_{hi}^2{w}_{hi}(x-x_{hi})+S\\ {\tau }_i\frac{d{a}_i}{dt}& =F_i(x)-a_i\forall i\in \{he,hi,le,li\}\end{array}\}$$

$$F_i(x)=\frac{1}{1+\exp (-{\Delta }_i(x-{\tilde{x}}_i))},$$

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_2}{dt}$$

$$H(x_1,x_2)=\frac{1}{2}(x_1^2+x_2^2),$$

$$\frac{d\mathbf{u}}{dt}$$

$$\displaystyle{\bf F}_{1}({\bf u})=\left(\begin{array}[]{c}\frac{-x-\sum_{i\in\{he,li,le\}}a_{i}w_{i}(x-x_{i})-a_{hi}^{2}w_{hi}(x-x_{hi})}{\tau_{x}},0,0,0\end{array}\right)^{T},$$

$$\displaystyle{\bf F}_{2}({\bf u})=\left(\begin{array}[]{c}0,\frac{F_{hi}(x)-a_{hi}}{\tau_{hi}},\frac{F_{le}(x)-a_{le}}{\tau_{le}},\frac{F_{li}(x)-a_{li}}{\tau_{li}}\end{array}\right),\;{\bf S}=\left(\begin{array}[]{c}\frac{S}{\tau_{x}},0,0,0\end{array}\right)^{T}.$$

$$f_1(x,\mathbf{a})$$

$${\mathbf{f}}_2(x,\mathbf{a})$$

$$\begin{array}{rl}{x}^{n+1}& =x^n+\Delta t{f}_1(x^n,{\mathbf{a}}^n)+\Delta t\frac{S}{{\tau }_x}\\ {\mathbf{a}}^{n+1}& ={\mathbf{a}}^n+\Delta t{\mathbf{f}}_2(x^{n+1},{\mathbf{a}}^{n+1})\end{array}\}$$

$$\begin{array}{rl}{x}^{n+1/2}& =x^n+\frac{\Delta t}{2}{f}_1(x^n,{\mathbf{a}}^n)+\frac{\Delta t}{2}S,\\ {\mathbf{a}}^{n+1}& ={\mathbf{a}}^n+\Delta t{\mathbf{f}}_2(x^{n+1/2},{\mathbf{a}}^{n+1}),\\ x^{n+1}& =x^{n+1/2}+\frac{\Delta t}{2}{f}_1(x^{n+1/2},{\mathbf{a}}^{n+1})+\frac{\Delta t}{2}\frac{S}{{\tau }_x}\end{array}\}$$

$$\begin{array}{rlr}& & \frac{d\mathbf{u}}{dt}=\mathbf{F}(\mathbf{u})\\ & & \mathbf{u}(0)={\mathbf{u}}^0\end{array}\}$$

$$\Vert \mathbf{u}(0)-\mathbf{u}(\tilde{t})\Vert =0$$

$${\mathbf{u}}^{n+1}={\mathbf{u}}^n+\frac{\Delta t}{2}\left(\mathbf{F}({\mathbf{u}}^{n+1})+\mathbf{F}({\mathbf{u}}^n)\right)$$

$${\tilde{\mathbf{u}}}^{n+1}={\mathbf{u}}^n+\frac{\Delta t}{2}\mathbf{F}({\mathbf{u}}^n),$$

$${\mathbf{u}}^{n+1}={\tilde{\mathbf{u}}}^{n+1}+\frac{\Delta t}{2}\mathbf{F}({\mathbf{u}}^{n+1}),$$

$${\mathbf{u}}_i^{n+1}={\tilde{\mathbf{u}}}^{n+1}+\frac{\Delta t}{2}\mathbf{F}({\mathbf{u}}_{i-1}^{n+1}),$$

$$\mathbf{F}(\mathbf{u},\mathbf{v}):={\mathbf{F}}_1(\mathbf{u})+{\mathbf{F}}_2(\mathbf{v})+\mathbf{S}$$

Full text

Keywords: Hodgkin-Huxley Type model, iterative solver methods

AMS subject classifications. 35K25, 35K20, 74S10, 70G65.

1 Introduction

We are motivated to model a nonlinear dynamic problem in neuroscience. The most prominent system to describe the dynamics of neural cells is the Hodgkin Huxley model [13]. It is characteristic for this class of models to exhibit highly nonlinear oscillations in response to some external input [2]. Sometimes we can observe chaotic oscillations, as for example in a small regime within the originally given parametrization of the Hodgkin and Huxley’s model [9]. Many subsequent models for biological oscillators have been either derived from this system or inspired by it. For details see [25] and [16].

To study such delicate nonlinear dynamics, it is important to deal with stiff ODE solvers, which preserve the structure of the solution, see [2] and [10]. Based on the high quality of explicit and implicit time-integrators, which can be combined with conservation scheme, see [29], we propose novel semi-implicit iterative methods, see [6].

The paper is outlined as follows. The model is introduced in Section 2. In Section 3, we discuss the different numerical methods and present the convergence analysis. The numerical experiments are done in Section 4 and the conclusion is presented in Section 5.

2 Mathematical Model

The classical Hodgkin-Huxley model is a parabolic partial differential equation with nonlinear reaction parts, see [13]. It models the dynamic behaviour of the the giant squid axon, which is a part of a neural cell. Neural cells transfer information with the help of voltage peaks (so called action-potentials). The voltage peaks base on the imbalance of the inner and outer ions and their diffusion, which is controlled by the potential difference across the cell’s membrane. The involved ions and channels are dependent on the type of neuron.

The standard Hodgkin-Huxley model is based on the flux of $Na^{+}$ and $K^{+}$ ions trough ion channels in the cell’s membrane and proton pumps to provide a non-equilibrium environment. Proton pumps move $Na^{+}$ ions out and the $K^{+}$ in by consuming ATP, forcing an imbalance of $Na^{+}$ and $K^{+}$ ions in the extracellular space (cell’s outside) and intracellular space (cell’s inside) respectively. The activation state of these ion channels is controlled by voltage (potential) at the membrane. When enough voltage is present, then the fast $Na^{+}$ channels start to open, launching a diffusion-driven inflow of $Na^{+}$ ions from the extracellular space outside into the inner-cell, changing the cell’s membrane potential towards positive values. After a short time the $Na^{+}$ channels close and keep this closed state over a short time (they are called to be refractory). Slow $K^{+}$ channels open delayed to the fast $Na^{+}$ channels, such that we have an outflow of $K^{+}$ ions, parallel to the closing $Na^{+}$ channels. This mechanism again introduces a change in the potential, turning it back to the initial potential. Putting these ideas together and taking into consideration, that the surface of neurons is geometrically rather complex and that the ion channels are not perfectly equal distributed over this surface, we can derive a partial differential equation to describe these spikes along an axon (see [13] or [16]):

[TABLE]

where $D_{m}\frac{\partial^{2}V}{\partial x^{2}}$ is the longitudinal conductivity. $I_{Na}$ are $I_{K}$ are the natrium and kalium induced currents, $I_{L}$ is the leak current and $I$ is some externally applied current. $n$ models the slow $K^{+}$ channel activation, while $m$ and $h$ describe the $Na^{+}$ channel activation and inactivation. All parameters can be determined experimentally.

Further we have

[TABLE]

and the transition rates:

[TABLE]

This system is original Hodgkin-Huxley PDE [13]. Now assuming an ideal model of a neuron (more specifically its axon) as a cable, such that the spatial sizes are homgeneous and independent, we can reduce the model to a system of ordinary differential equations of the form (see also [13]):

[TABLE]

There also exist model-reductions of the HH model, mostly based on 2D ODEs (e.g. see [17, 25]). One of the most famous one is the FitzHugh-Nagumo (FHN) model [4, 23]. Such models cannot show chaotic behaviour as a consequence of the Poincare-Bendixson theorem [12]. The FHN model can also be interpreted as a generalisation of the Van-der-Pol Systems and is given as:

[TABLE]

2.1 Hodgkin-Huxley Type Models

To the best of our knowledge there exists no formal definition of which models exactly belong the class of Hodgkin-Huxley type systems. Informally we refer to Hodgkin-Huxley type systems as differential equations as a special class of potentially nonlinear oscillating systems, where oscillations of an observable quantity are induced by the interplay with some independent but dynamic activations. Closest to a definition of this class is the generalized deterministic Hodgkin-Huxley equation by Tim Austin [1]. Based on the definition given from [1] and observations we propose the following definition for the class of Hodgkin-Huxley type (HHT) systems

[TABLE]

which can be interpreted as a special case of reaction-diffusion systems.

Here $o$ describes an observable quantity, $I$ describes some external input function and $a_{i}$ are activation quantities. $f_{o}$ couples the observed quantity to the activations and may contain partial differential and integral operators, while the $f_{i}$’s describe analytic couplings of the activation back with the observable quantity. From a modeling perspective we can sometimes an ideal case, where the spatial domain of PDE (1) is homgeneous and independent, such that the system reduces to the following ODE:

[TABLE]

This way our system contains naturally well-known ODEs used to model neural dynamics. To the best of our knowledge this preliminary definition contains most systems which has been attributed as Hodgkin-Huxley typed in scientific literature so far. We are aware that a special class of systems is not directly captured, the hybrid dynamical system models, since their solutions are discontinuous (to speed up computations of trajectories) [18], although the underlying continuous part of the system is.

2.2 A Hodgkin-Huxley Type Nonlinear Disease Dynamics Model

Trough this paper we deal with a HHT model appearing in neural modeling from neuroscience [14] and as the deterministic part for a stochastic disease model in neuropsychiatry [15]. We chose this model as a representative system for the class of HHT models because it exhibits a rich amount behaviour in response to a constant input. The model is given by the following system of ordinary differential equations

[TABLE]

clearly fitting in the HHT class defined in the equations (3-4). We have x as an observable, where peaks represent events within the disease. Further $\{he,hi,le,li\}$ are the different activation types, operating on two time scales. Elements starting with $h$ describe the fast time scale and model a high activation threshold, while elements starting with $l$ describe the slow time scale with low activation threshold respectively. $e$ describes an excitatory and $i$ a corresponding inhibitory quantity. $F_{i}$ are sigmoidal functions of the form

[TABLE]

where $\tilde{x}_{i}$ is the half-activation levels and $\Delta_{i}$ is the steepness of the sigmoidal function. The fast excitatory quantity is assumed to activate instantaneously, so the model always has $\tau_{he}=0$, implying $a_{he}=F_{he}(x)$. As a consequence we also reduced the dimension of our dynamical system from 5 to 4.

3 Numerical Methods

In the following we apply and discuss composition methods, as well as structure preserving methods based on finite difference and iterative schemes, which are also known to be successful in approximating solutions for various reaction-diffusion type equations. We restrict us to composition methods, while also in the literature, there exists different other types of solver methods, e.g., tailored multi-step methods as the Rush-Larson method, see [24].

For notational simplicity we assume S to be time-independent such that the system gets autonomous. Further we introduce the following notation:

• •

${\bf u}=(x,a_{hi},a_{le},a_{li})^{T}=(x,{\bf a})^{T}$ is the exact solution,

• •

${\bf u}^{n}=\left(x(t^{n}),a_{hi}(t^{n}),a_{le}(t^{n}),a_{li}(t^{n})\right)^{T}=(x(t^{n}),{\bf a}^{n})^{T}$ is defined as the solution at the time-point $t^{n}$,

• •

Analogously ${\bf u}_{i}^{n}=(x_{i}(t^{n}),{\bf a}^{n}_{i})^{T}$ is defined as the iterative solution of ${\bf u}$ in the $i$-th iterative step at the time-point $t^{n}$.

Bold letters indicate vectorial objects and $(\;\cdot\;)^{T}$ is the transpose.

3.1 Composition with respect to Hamiltonian Systems

If we apply a Van-der-Pol oscillator, which is a very simple Hodgkin-Huxley type system, we can reformulate the oscillator with respect to the non-stiff case into a Hamiltonian system and apply splitting approaches for the Hamiltonian systems. The Van-der-Pol oscillator is given as:

[TABLE]

where for $\mu=0$, we obtain the harmonic oscillator with the Hamiltonian system

[TABLE]

although also other approaches are possible to uncover the systems hamiltonian [26]. With these structural observations the idea is to apply such composition methods, which are known for the Hamiltonian system, i.e. Semi-implicit Euler scheme and Störmer-Verlet scheme [10, 11]), which are symplectic schemes if they are applied to a Hamiltonian system.

We introduce the following composition in operator notation for the disease model:

[TABLE]

where

[TABLE]

Basing on this we define

[TABLE]

such that the algorithms are given as:

• •

Semi-implicit Euler scheme:

[TABLE]

• •

Störmer-Verlet scheme:

[TABLE]

Remark 1

We can solve equations depending explicit on an $x$ and implicit on ${\bf a}$ directly, since the equations can be trivially rearranged on account of the linearity and independence on ${\bf a}$ in ${\bf f}_{2}$.

Remark 2

For the semi-implicit Euler we have a global convergence order of $\mathcal{O}(\Delta t)$ and for the Störmer-Verlet $\mathcal{O}(\Delta t^{2})$.

3.2 Iterative Schemes Based on Finite Difference Schemes

We deal with the disease model, which is given as:

[TABLE]

We assume to deal with a system containing exactly one periodic orbit (in properly parameterized regime). This implies there exists a $\tilde{t}>0$ such that for all points ${\bf u}_{0}$ starting on this orbit holds:

[TABLE]

We call the smallest $\tilde{t}$ the period of an orbit. We apply a semi-impicit Crank-Nicolson scheme (CN), see also [29], which is conservative and given as:

[TABLE]

Here, we have a nonlinear equation system, which have to apply additional nonlinear solvers, e.g. Newton’s method. Therefore, we propose iterative schemes, which embed via iterative step to the semi-implicit structures.

Remark 3

The semi-implicit CN method can be derived via operator-splitting approach:

[TABLE]

where the first equation (3) is explicit and can be done directly, the second one (3) is implicit and solved with a fixpoint scheme as:

[TABLE]

where the starting condition is ${\bf u}_{0}^{n+1}={\bf u}^{n}$ and we apply $i=1,\ldots,I$, while $I$ is an integer and we stop if we have the error bound $\left\lVert{\bf u}_{i}^{n+1}-{\bf u}_{i-1}^{n+1}\right\rVert\leq\varepsilon$ with $\varepsilon$ as an error bound.

3.2.1 Semi-implicit Integrators

In the following, we deal with semi-implicit integrators. We introduce the following the following convention for intermediate results:

• •

We initialize the iterative scheme with the solution in time point $t^{n}$, i.e. ${\bf u}_{0}^{n+1}={\bf u}^{n}$.

• •

We set the approximation for the next time point $t^{n+1}$ with the iterative solution in the $i$-th iterative step, i.e. ${\bf u}^{n+1}={\bf u}_{i}^{n+1}$

• •

We will denote the splitting from equation (6) as follows:

[TABLE]

We compute the approximations ${\bf u}(t^{n})$ at the time points $n=1,2,3,\ldots,N$ coupled with a fixed-point iteration, where $t^{N}=T$. The initialization of the iterative scheme is given with the initial condition of the equations (11) as $u^{0,1}=u^{0}$. For now the time step is defined as $\Delta t:=t^{n}-t^{n-1}$, while the error bound is given as $\varepsilon$. Based on this information we define the first three solvers with algorithms (1-3).

Remark 4

The semi-implicit CN scheme based on the iterative approach is asymptotical conservative [7].

Further we define two multipredictor multicorrector methods with algorithms (4) and (5).

3.3 Adaptive Time Step Control of the Iterative CN Scheme

To improve the numerical results in the critical time-scales (i.e. the stiff parts of the evolution equation) we apply adaptive time step approaches. We define the following norms:

• •

Absolute norm:

[TABLE]

• •

Maximum-norm:

[TABLE]

The relative error is given as:

[TABLE]

3.3.1 PID-Controller

We apply the following simple error-estimate (see [21]), where we compute the time step for a given tolerance $\varepsilon$ at a timepoint $t^{n}$:

[TABLE]

where we assume the emprical PID (Proportional-Integral-Differential) parameters $k_{P}=0.075,k_{I}=0.175,\;k_{D}=0.01$. For the initialisation, means for $n=1$, we only apply the $I$ part, while for $n=2$ we apply the $I$ and $P$ part and for all later time steps (where we have all the parts $e(t^{n-2}),e(t^{n-1}),e(t^{n-2})$), we apply $I,P,D$.

3.3.2 Classical Time Step Controller for the ICN

We apply an additional automatic time step control which is given as following with a two scale ansatz, where we compute an approximation via large step $\Delta t$ and compare the solution with m consecutive substeps of length $\frac{\Delta t}{m}$ to give another approximation, which should be close to the large step if the approximator is accurate enough, given the current time step. Solutions are rejected until the time step is small enough, which implies the approximation error is smaller than some bound. We apply the following time step controller for second order schemes:

[TABLE]

where $\Delta t^{*}$ is the optimal time step while ${\bf u}_{\Delta t}$ is the approximation by applying m small time steps an ${\bf u}_{m\Delta t}$ is the solution of an equivalent length large time step.

3.4 Time Step Controller for the Runge-Kutta Methods

We extend the multipredictor-multicorrector algorithm of order $4$, see Algorithm (4) and an iterative CN+Simpson-Rule of order $4$, see Algorithm (5):

Lemma 1

We deal with 4th order time-integrator methods with tolerance $\varepsilon$. Further, we assume that we have a 4th order numerical solver, which is give as ${\bf u}(t+\Delta t)=A_{\Delta t}\;{\bf u}(t)$ and ${\bf u}(t)$ is the exact solution at time $t$. We apply the $||\cdot||_{p}$-norm as a given vector norm, e.g., in the Banach-space.

Then the adaptive time stepping is given as:

[TABLE]

Proof

We assume $\left\lVert{\bf u}-{\bf u}_{\Delta t}\right\rVert=\varepsilon$, which is a prescribed tolerance.

We apply 2 different time-steps:

• •

A single large time-step $\Delta t$ with:

[TABLE]

• •

A multiple small time-step $\Delta t/m$ with:

[TABLE]

The local truncation error is given as:

[TABLE]

and we assume to have the approximation:

[TABLE]

which can be interpreted as a scaling of the error estimates.

Using the norm property we can now pull out the divisors:

[TABLE]

we can simplify the divisors:

[TABLE]

we assumed \left\lVert\mbox{\boldmathu}(t^{n})-\mbox{\boldmathu}_{\Delta t^{*}}(t^{n})\right\rVert_{2}=\varepsilon, which is our error control, such that we obtain the following crude approximation:

[TABLE]

Then, the adaptive time stepping is given as:

[TABLE]

The improved automatic time step controlled $4$-th order methods are now given with algorithms (8) and (9).

4 Numerical Results

Trough this section we present a short analysis of the dynamical system in combination with the performance of the in previous section derived solvers. For the implementation we used Julia111https://julialang.org/ 1.1. A Jupyter notebook containing the implementation of this section can be found online under https://git.noc.ruhr-uni-bochum.de/ogierdst/solving-hodgkin-huxley-type-systems/.

We deal with the disease dynamics model (2) and the parametrization taken from [15]:

[TABLE]

Note that since $\tau_{he}=0$ we obtain a reduced system of order 4, where $a_{he}=F_{he}(x)$. This choice corresponds to an instantaneous activation of $a_{he}$, effectively reducing the system’s dimension to 4.

4.1 Exploring Structural Properties via Computational Bifurcation Analysis

We start by exploring the system’s overall behavior for varying $S\in[0,400]$. This section is not ment to replace a rigorous dynamical system analysis, but to outline its coarse structure to ease the analysis of the solvers. For convenience we use Tsit5 from the JuliaDiffEq package [19] as the solver when not otherwise stated. This way we provide a tested baseline as a foundation to compare the implementation of our solvers to.

As a first step we extract the system’s fixed-points, which are given by setting the change in all dimensions to zero. Formally we first rewrite the model (2)

[TABLE]

and set it to zero, i.e.

[TABLE]

Here ${\bf u^{*}}$ denotes a fixed point. The system’s special structure allows us to reduce this problem to one dimension, as

[TABLE]

which results in

[TABLE]

It can be easily shown that this function is unbounded and strictly monotonically decreasing for our chosen parametrization. This implies that there is a single fixed point for each S. We obtain the corresponding $a_{i}^{*}$’s explicitly by plugging the solution back into the corresponding equations. Approximating some fixed points with Newton-Raphson and linearising around these gives an idea of its stability properties. This yields the Jacobian $J_{ij}=\frac{\partial f_{i}}{\partial u_{j}}|_{{\bf u}={\bf u}^{*}}$, which is explicitly:

[TABLE]

Note that $f_{i}$ is the disease models i-th equation while $F_{i}$ denotes the sigmoidal function for the corresponding activation.

Further we approximate the Lyapunov spectrum as a measure for the divergence of nearby trajectories to obtain information about the system’s stability properties. The Lyapunov spectrum is formally defined as

[TABLE]

where $\alpha_{i}$ are the eigenvalues of $M(t)M^{T}(t)$. Here $M$ denotes the discrete time evolution operator. We carry out the numerical approximation of the lyapunov spectrum with ChaosTools [3]. The results are presented in figure 2. Lyapunov exponents can be seen as a simple characterization for the stability of manifolds, where a negative exponents indicate attraction, positive exponents indication repulsion and an exponent of zero indicates conservation.

These figures together suggest a Andronov-Hopf bifurcation around $S\approx 20$, where in the interval $[0,~{}20)$ the fixed point is a stable one. After this we see a maximal Lyapunov coefficient of value zero paired all other coefficients negative, which is associated with stable cycling. Around $S\approx 180$ we see that the largest Lyapunov coefficient gets positive. This is possibly associated with the onset of chaos. Further around $S\approx 340$ the systems gains stability again, which is in turn possibly associated with the end of chaotic behavior, returning to stable cycling again. Around $S\approx 100$ we see the two real eigenvalues becoming complex. We failed to associate this observation with any phenomenon.

Now that we have worked out the coarse system structure we move on to confirm details computationally. We start by approximating solutions for arbitrary S from each identified interval, namely $[5,100,180,255,340,400]$, with algorithm 2 with tolerance $\varepsilon=10^{-7}$, time step $\Delta t=0.01$ and the maximum number of iterations $I=10$. The results are visualized in figure (3). It can be clearly seen that for $S=5$ the fixed point is attracting, while all other choices of S yield oscillations, which is on par with the previous computational analysis of the Jacobian and Lyapunov spectrum. The choice $S=255$ suggests either an unstable solver or chaotic cycling behavior. Please note also that solver takes up some time to settle, i.e. moving from the initial condition into an orbit.

With this basic structural guesses we move forward towards a computational bifurcation analysis, as to the best of our knowledge no analytic work is available about the general structural properties of Hodgkin-Huxley type systems and especially our disease dynamics model. We will use two related techniques to quantify the systems behavior computationally, namely Poincaré maps and interspike intervall (ISI) distributions. For both techniques we will use the same section. This will also give us some clues about very basic stability and correctness properties of the in the previous section constructed solvers.

Poincaré sections allow us to study the behaviour of continuous high-dimensional system with a geometric description in a lower-dimensional space, see [28]. The basic idea is to reduce the system to a continuous mapping $T$ of the applied plane $S$ into itself, means we have:

[TABLE]

Therefore we reduce the continuous flow into a discrete-time mapping. The Poincaré section of the hyperplane $<(1,0,0,0),u>=40$ can be seen in figure (4).

A closer look into the regions with the first branch and the last merge reveals can be found in figure (5). The found structures can be identified as classical period doubling and period halving, which are more pointers towards the existence of chaotic behavior, as they usually indicate the onset and the end of chaotic regimes. Computing the position of three branching points trough a finer step size for S in the Poincaré section, starting with the second branching point (i.e. $\approx(175.1,178.8,179.6)$, yields a ratio close to Feigenbaum’s constant, suggesting period doubling. The same structure can be found on the other side at the end of the hypothetically chaotic regime, suggesting period halving (i.e. $\approx(342.25,340.0,339.5)$). While a rigorous analysis is out of the scope of this paper, we take the worked out arguments to support the assumption, that chaotic cycling is actually present as a property of the dynamical system and not as a numerical artifact of instabilities in the used solvers.

As a next step we generate the interspike interval distributions for the same section, which is basically the distribution of time between two consecutive intersections of this plane of the solution, which starts in the corresponding attractor. This distribution is approximated by fixing S and solving the system for a fixed time interval (here [0,10000]). The bifurcation plot for each in this paper defined scheme can be found in figure (6).

4.2 Convergence Study

Now we test the convergence behaviour of the schemes with fixed time step. We arbitrarily take one configuration of S for the stable as well as the chaotic cycling, namely $S\in\{100,253\}$. The convergence analysis is conducted as follows. We start the first approximation with initial condition $u(0)=(0,0,0,0)$ and a time step $\Delta t=0.5$. With each consecutive approximation we reset the inditial condition and halve the time step while fixing all other parameters. The fixed parameters are $\varepsilon=10^{-7},I=5$. Consecutive approximations are now compared at the overlapping time points by integrating over the difference of these consecutive approximations. The results can be found in figure (7).

If the error shrinks with each halving the scheme converges to a solution, which should in the case of stable cycling be the corresponding solution of our system. In the case of chaotic cycling the solution converges only for this specifically given time interval, as in the presence of chaos nearby trajectories diverge with exponential speed. This divergence cannot be handled in general by our solvers for long time scales. On the one hand we usually cannot hit a solution exactly with only machine precision available, which may already be another solution trajectory which diverges exponentially. On the other hand we can, again by machine precision limited, not reduce the time step for an arbitrarily large time interval, as computations get unstable for too small time steps.

4.3 Analysis of the Adaptive Schemes

Finally we want to evaluate the performance of the adaptive schemes, namely algorithms 6, 7, 8 and 9, side by side with the fixed time step schemes. For this we start by evaluating the following local norms:

• •

$L_{2}$-norm of the solutions for each time-points:

[TABLE]

• •

$L_{2}$-norm of the derivations for each time-points:

[TABLE]

which result in the following global norms:

• •

$L_{2}$-norm of the solutions in the time domain:

[TABLE]

• •

$L_{2}$-norm of the derivations for each time-points:

[TABLE]

where the derivations are approximated as $\frac{dx(t^{n})}{dt}\approx\frac{x(t^{n+1})-x(t^{n})}{\Delta t}$. Again we chose two S arbitrarily (here {100, 258.15}) from the stable and chaotic cycling regions. The solvers are configured as follows:

[TABLE]

For the PID-controller we hand-tuned the parameters to the following values:

[TABLE]

. The time domain is $[0,10000]$. The local norms can be found in figure (8) while the global norms are listed in table (2) side by side with benchmarked runtimes. For the stable cycling we observe that all solvers nearly agree on the given time interval. Only the ISIE (1) and ISV (3) schemes are a bit off. No agreement is found in the chaotic case.

We further capture statistical features of the adaptive schemes fluctuations by computing the expectation and variance as follows, assuming that the same point does not lie on the approximation twice for our chosen time intervals:

• •

Expectation:

[TABLE]

• •

Variance:

[TABLE]

where $t^{n_{1}}<t^{n_{2}}$. The results can be found in table (3). We can take from these tables that the adaptive RK4 schemes (8 & 9) can handle larger time steps while bounding the local error. Note that this does not help in the assumed chaotic case, as in chaos nearby trajectories diverge with exponential speed. Still, with a small enough error bound we are able to somewhat bound the global error for small time intervals.

5 Conclusion

We started by giving a definition for Hodgkin-Huxley type systems and some characteristics. Informally we refer to Hodgkin-Huxley type systems as differential equations of a special class of potentially oscillating reaction-diffusion type systems with local activation and inactivation mechanisms. Based on different assumptions about these systems we derived some solvers and took an characteristic example, whose structure has been analyzed computationally. We found a potential period doubling and period halving around an unstable regime, which we assume to be chaotic. This chaotic regime has been examined further computationally, exposing a spiral structure via special Poincaré section from computational neuroscience called interspike interval bifurcation, which is not visible in the vanilla Poincaré section. The hereby taken approach can be seen as a basic framework for to guide numerical analyses of solvers for Hodgkin-Huxley type systems.

All solvers agree on the basic spiral structure, excepting the PIDICN which was unable to unravel higher windings, yielding much noise across the diagram in this area for the given parametrization. We also observed that the computations of the interspike interval bifurcations with adaptive higher order solvers lead to less noisy looking structures. This gives us pointers that these solvers, while not agreeing on solutions due to the potential chaos, still somewhat preserve the character of solutions.

In the future we look forward to analyze stochastic definitions of Hodkin-Huxley type systems. We also plan to examine geometrical and dynamical properties of the interspike interval bifurcation rigorously, providing a better foundation to understand the properties of numerical solvers for these kind of systems.