Skip to content

TBMT42 Lab 1B: Analysis of non-linear properties of a system

This laboratory exercise will introduce a couple of non-linear properties and different non-linear systems. Firstly, we will analyze transient model behaviours vs stationary model behaviours and how they can change over time. Secondly, we will investigate the different types of equilibriums there are and how domains of attraction influence our model's ability to find these equilibriums. Lastly, we will apply the knowledge from the prior parts to perform a bifurcation analysis. The bifurcation analysis is used to study qualitative behaviours of models and is very useful when working with nonlinear systems, which might have unforeseen (or hard to predict) behaviours built into them.

Throughout the lab you will encounter different "admonitions", which as a type of text box. A short legend of the different kinds of admonitions used can be found in the box below:

Different admonitions used (click me)

Background and Introduction

Useful information

Guiding for implementation

Here you are tasked to do something

Reflective questions

This module will give you pre-made scripts for Python. You can download the example files here. You are free to use any software, however, if you don't use Python you will need to redo the pre-made scripts yourself. For installation and setup of, see Python installation and Python packages. The packages required for this lab is found in requirements.txt. Note that the pre-made scripts are using the SUND package.

Introduction

First thing to do, is to read the introduction "hidden" below.

Non-linear systems

One fundamental aspect of biological systems is their inherent complexity. Biological processes often involve intricate interactions between multiple components, feedback loops, and non-intuitive dynamics. In many cases, linear models fail to capture the full complexity and richness of these systems. This limitation has led to the recognition of the importance of nonlinear models in understanding and describing biological phenomena.

Nonlinear systems play a crucial role in biology, as they accurately represent the intricate relationships and interdependencies among different variables. Unlike linear systems, where the output is directly proportional to the input, nonlinear systems exhibit more intricate behaviours, often displaying unexpected patterns and emergent properties. Nonlinearity arises when the relationship between variables is nonlinear, meaning that small changes in the input can result in disproportionate and non-proportional changes in the output. In biological systems, nonlinear dynamics can arise from various factors. For example, the interactions between biomolecules within a cell, the regulatory networks governing gene expression, or the feedback mechanisms in neuronal networks all give rise to nonlinear behaviours. Nonlinear systems can exhibit phenomena such as oscillations, chaos, bifurcations, and stability transitions, which are critical for understanding the behaviours of biological systems at different scales.

To study nonlinear biological systems, mathematical models need to incorporate nonlinear relationships between variables. These models often utilize differential equations, difference equations, or stochastic equations to describe the dynamic behaviours of the system. Nonlinear models can capture a wide range of phenomena, such as the growth of populations, the spread of infectious diseases, the dynamics of biochemical reactions, and the behaviour of neural networks.

This laboratory exercise will introduce a couple of non-linear properties and different non-linear systems.

Transient stationary behaviours

Concept of transient and stationary behaviours

Firstly, let's delve into the concept of transient and stationary behaviour. These two phases describe the different states a system can exhibit. The system can undergo adjustments, such as responding to an external stimulus, which is known as transient behaviour. On the other hand, the system can be in a stationary or steady-state phase where it remains unchanged over time. Since biological systems are often subjected to external stimuli, they frequently transition between these phases. The duration of this transition, known as the transient time, may vary between systems, resulting in different amounts of time spent in the stationary phase.

In this part of the lab, we will work with a simple model describing the signalling between glucose and insulin. As this system is constantly active, and can be affected by external factors, it is suitable to observe transient and stationary behaviours.

Glucose-insulin model

Consider a mathematical model that describes the interaction between glucose and insulin in the human body. 

d/dt(G) = k1 * Gb     - k2 * I 
d/dt(I) = k3 * Ib * G - k4 * I

Here, \(k1\) represents the basal rate of glucose production, and \(k2\) represents the insulin sensitivity or the efficiency of insulin in promoting glucose utilization. Additionally, \(k3\) represents the rate of insulin secretion in response to glucose levels, and k4 represents the rate of insulin degradation. \(Gb\) and \(Ib\) are the basal levels of glucose and insulin.

You have been given the model files in the example files, you should now make sure that you can run the prepared files and find the stationary phase for the glucose-insulin model.

1. Stationary behaviour

Now, find the stationary phase for the model file you have been given, see Models/glucoseInsulin.txt. You can simply do this by running the TransientStationarybehaviour.py script you have been given.

Have we reached a stationary phase?

If the derivatives of the states equal zero, i.e. the state values are non-changing, then we have reached an equilibrium.

Moving on, you will try to move the model out of its stationary phase into a transient phase.

Transient behaviour

Now, let us imagine a scenario where an individual consumes a high-carbohydrate meal. This leads to an increase in blood glucose levels, triggering the release of insulin from the pancreas. Initially, the insulin response may be delayed, causing a temporary mismatch between the rising glucose levels and insulin secretion. As a result, the blood glucose levels continue to rise before they start to decline. We can add a meal to this system by introducing a meal_influx, which adds glucose to the system over a period.

To do so, we must introduce an external effect that can affect the model - this could be the uptake of glucose from a meal.

The meal effect implementation

As you might have observed, the meal effect is already present as the meal_influx input in the model file. This is one way of implementing time dependent effects, where a constant can be used to set the value of the input effect from the simulation call. This allows us to simulate with and without the meal effect using meal_influx. This could be done using a parameter as well, but you would need to exclude this parameter from the optimisation call. Below you see the state implementation.

########## STATES
d/dt(G) = k1 * Gb     - k2 * I + meal_influx
d/dt(I) = k3 * Ib * G - k4 * I

########## INPUTS
meal_influx = meal_influx @ 0

Here, meal_influx is set to be 0 per default. We can create activity objects to govern how this constant should be changed over time. 

activityMeal = sund.Activity(timeunit='m')
activityMeal.AddOutput(sund.PIECEWISE_CONSTANT, "meal_influx", tvalues = [-100, 0, 120], fvalues=[0, 0, 0.1, 0], feature=True)

where, meal_influx=0.1 if 0<t<120, else meal_influx=0

This has been implemented in the TransientStationaryBehaviour.py for you.

Now we want to find if the model can find a new stationary phase with this meal effect.

2. Find new stationary phase

You should now experiment with different meal-durations for the simulated meal. Start by uncomment the code segment Meal effect. The script is set up for you so you can simply change the timeDuringMeal variable in the TransientStationaryBehaviour function to construct the time vector.

Have we reached a stationary phase?

Let us look at the 2nd Stationary phase figure. Again, if the derivatives of the states equal zero, i.e. the state values are non-changing, then we have reached an equilibrium. Since we have added a continuos input, the models should stabilize at a new equilibrium point, as the insulin secretion will compensate for the added influx of the glucose.

If we observe the model’s behaviour, we can see that the changes applied to the model are evened out and a new stationary equilibrium is found. The model enters a transient phase to move between the two stationary phases. This should also be true when the effect of the meal disappears, and the model no longer has an external effect.

3. Return to the original stationary phase

To return to original stationary phase we need to simulate further, after the meal has ended, and we do this by setting the timAfterMeal variable to be greater than timeDuringMeal variable, also remember to uncomment the Post meal section. How long does the model need to return to the original baseline? Are the timings the same for moving from the first to the second equilibrium and vice versa (compare the 2nd Stationary phase and Return to original stationary phase figures)?

Transient duration

Is the transient duration, or transient time, the same for moving between the first and second equilibriums?

You have now explored and found two different equilibrium points and you should be able to find more if you alter the impact of the glucose influx from the meal. However, are these equilibriums the same? The short answer is no, they are not. Equilibrium's can be stable or unstable, which can have big impacts on how we use and work with nonlinear systems. In the section below you can read a bit more in the difference between stable and unstable equilibriums.

Stable vs unstable stationary phases

With the extra influx of glucose due to the meal, the system reaches a new equilibrium or steady-state condition. In this stationary phase, the body's insulin secretion and glucose uptake mechanisms adjust to the high glucose levels. Insulin helps facilitate the uptake of glucose by cells, leading to a reduction in blood glucose levels. Eventually, the insulin secretion and glucose uptake reach a balance, resulting in a stable blood glucose level within a certain range. Continuing, the meal_influx ended, and the system will return to the same balance it had prior to the meal.

This showcases stable vs unstable steady states. The glucose-insulin dynamics were stable for a period when the meal was having an effect. However, as the meal ended the system quickly left this equilibrium. This is due to the system being in an unstable equilibrium point, where the system per definition is in balance but is highly subjective to external changes. We can illustrate this using the figure below. Here, the ball is in a stationary phase on top of the hill. However, it can easily be offset and roll down the hill. This is the transient phase. Ultimately, it ends up at the bottom of the hill. This on the other hand, is a stable stationary phase where it will return even though external factors may act upon the ball.

Stable vs unstable stationary phases

Domains of attraction

To be able to properly work with stable and unstable equilibriums and understand how the model gravitates towards them, we need to understand if they are attracted to the domain of the equilibrium.

Domains of attraction

A domain of attraction is associated with a particular equilibrium point or steady state of the dynamical system. These points represent a stable configuration of the model states, the derivatives equal 0, that the system tends to approach over time. The domain of attraction of an equilibrium point comprises all the initial conditions from which the system's trajectory converges towards that equilibrium, meaning all starting positions that will end up in the equilibrium.

The shape and size of the domain of attraction depend on the stability properties of the equilibrium point. Stability is often classified into three types: stable equilibriums, unstable equilibrium, and semi-stable equilibriums. In the figure below, you can see a graphical illustration for a system with two parameters (\(\theta_1\) and \(\theta_2\)). The grey area is the domain of attraction, the orange circle represents the equilibrium, the blue square represents an initial condition, and the blue arrow represents the movement from that initial condition.

Domains of attraction

  • Stable equilibrium: If a small perturbation in the initial condition leads the system to return to the equilibrium point, it is called a stable equilibrium. The domain of attraction for a stable equilibrium is typically a region surrounding the equilibrium point.

  • Unstable equilibrium: If a small perturbation in the initial condition causes the system to move away from the equilibrium point, it is referred to as an unstable equilibrium. The domain of attraction for an unstable equilibrium may be empty or include only specific initial conditions that lead to transient behaviour.

  • Semi-stable equilibrium: In some cases, the stability of an equilibrium point can depend on the direction of perturbation. If the system returns to the equilibrium point when perturbed in one direction but moves away from it when perturbed in another direction, it is called a semi-stable equilibrium. The domain of attraction for a semi-stable equilibrium can be asymmetric and may differ based on the direction of perturbation.

A domain of attraction can provide valuable insights of the behaviour of biological systems. They can reveal information about the stability, resilience, and robustness of biological processes and help in predicting the long-term outcomes of dynamic models. By studying the domains of attraction, researchers can gain a better understanding of the dynamics of biological systems and make predictions about their behaviour under different conditions.

To start off with, we will observe the equilibrium for a model with a semi-stable equilibrium point, the SEIR model.

The SEIR model

The SIR model is commonly used to study the spread of infectious diseases. The SIR model divides the population into three compartments: susceptible (S), infected (I), and recovered (R). SEIR is an alteration of the SIR model, which also includes a compartment for exposed (E) individuals. This model framework was widely used for trying to predict the covid-19 spread and impact on the society, see examples i) SEIR model for COVID-19 dynamics incorporating the environment and social distancing, ii) An Enhanced SEIR Model for Prediction of COVID-19 with Vaccination Effect, and iii) SEIR modelling of the COVID-19 and its dynamics.

The equations for our SEIR model are as follows:

########## STATES
d/dt(S) = -B*S*f(I)/N
d/dt(E) =  B*S*f(I)/N - C*E/N 
d/dt(I) =  C*E/N - y*I
d/dt(R) =  y*I

########## FUNCTIONS
f(X) = max(X-5, 0)

Here, \(B\) represents the exposure rate, \(C\) represents the transmission rate, \(y\) represents the recovery rate, and \(N\) is population size. The function \(f()\) is a threshold function that controls that a disease outbreak cannot occur if there is less than 6 people infected.

4. Find the equilibriums

Run the script SemiStableEquilibrium. The script will plot the transition to a equilibrium point for the SEIR model, for two different initial guesses. Your job here is to observe the model structure and try to find an initial starting point that leads to a second equilibrium point.

Changing initial guesses

Do this by uncommenting the IC3 section and setting the values for IC3. Things to note:

  • The model assumes a total population of 1000 (N in the model)
  • We disregard the value of \(R\) in the equilibrium point, as all changes in transitioning to the equilibrium point will end up here.
Have you found a new equilibrium point?

In line with the previous section, you should have found a new equilibrium point where the derivatives equal 0.

  • What is the model behaviour that gives rise to these different equilibriums?
  • What is the interpretation of these equilibrium points on the population?
In what direction is the equilibrium stable?

As you have read, a semi-stable equilibrium point is stable in moving in one direction of the parameter space, while being unstable when moving in another direction. Following are some questions regarding this behaviour:

  • Which initial condition is crucial for the model to change between the equilibrium?
  • In what direction (the values of the initial condition) are the model able to stay at the semi-stable equilibrium?
  • In what direction will the model leave the semi-stable equilibrium and return to the stable equilibrium?

To put this into contrast, we will also evaluate the equilibrium points for the FitzHugh-Nagumo model.

The FitzHugh-Nagumo model

The FitzHugh-Nagumo model, which is often used to describe the electrical activity of excitable cells, such as cardiac cells and neurons. The FitzHugh-Nagumo model is a simplified version of the Hodgkin-Huxley model and consists of two coupled differential equations:

########## STATES
dV/dt = V - (V^3/3) - w + I
dw/dt = e * (V + a - bw)
dI/dt = 0

In these equations, \(V\) represents the membrane potential, \(t\) represents time, \(w\) represents a recovery variable, \(I\) is an applied current, and \(e\), \(a\), and \(b\) are parameters that control the dynamics of the system.

5. Finding an unstable equilibrium

The FitzHugh-Nagumo model is used to describe the electrical activity in excitable cells and exhibits an unstable equilibrium. Explore the model behaviour by testing different initial conditions in the file UnstableEquilibrium.

Have we found an unstable equilibrium?

An unstable equilibrium point is associated with needing a very specific initial guess(es) to reach it. Unlike earlier examples, a dynamic model could be in equilibrium while having a transient behaviour if that behaviour has a constant period time (see figure below). What happens with the model behaviour as we approach the unstable equilibrium?

Period time

Bifurcation analysis

In the last section of this lab, we will go through bifurcations and how they impact our systems. Read more about bifurcations below.

Bifurcations

Bifurcations refer to qualitative changes in the behaviour of a system as a parameter or set of parameters is varied. Bifurcations can lead to the emergence of new steady states, periodic oscillations, or chaotic behaviour in the system. In mathematical models, bifurcations occur when the equations governing the system undergo a qualitative change in their solutions due to changes in the model's parameters. This can result in a change in the stability, number, or types of equilibrium points, limit cycles, or chaotic behaviour exhibited by the system.

Bifurcations are important in biological modelling because they can help explain the observed transitions and changes in behaviour that occur in biological systems. They provide insights into how small changes in parameters can lead to significant shifts in the dynamics of the system.

There are several types of bifurcations that can occur in mathematical models. Some common examples include:

Saddle-node bifurcation

In this bifurcation, an equilibrium point collides with and disappears from the system as a parameter is varied. It leads to the creation or destruction of stable equilibrium points and can be illustrated using the equation \(\frac{dx}{dt} = r + x^2\). As \(r\) varies we get the following behaviors (illustrated in the figure below).

  • r<0: we have a stable equilibrium point (at \(-\sqrt{-r}\)) and an unstable equilibrium (at \(+\sqrt{-r}\)).
  • r=0: we have a stable equilibrium point.
  • r>0: we have no equilibrium points.

Saddle-node bifurcation

Pitchfork bifurcation

In this bifurcation, an equilibrium point undergoes a symmetry-breaking transition, resulting in the emergence or disappearance of two equilibrium points. In the first case, two stable equilibrium points emerges after the bifurcation point - this type is called a supercritical pitchfork bifurcation. In the other case, two unstable equilibrium points exist before the bifurcation and disappear when moving past the bifurcation point - this type is called a subcritical pitchfork bifurcation.

The showcase this, we will observe a simple equation that undergoes a supercritical pitchfork bifurcation - \(\frac{dx}{dt} = rx - x^3\). As \(r\) varies we get the following behaviors (illustrated in the figure below).

  • r<0: we have a stable equilibrium point.
  • r=0: we have a stable equilibrium point.
  • r>0: we have an unstable equilibrium point (at \(0\)) and two stable equilibrium points (at \(\pm \sqrt{r}\)).

Pitchfork bifurcation

Note that the unstable equilibriums only exist for \(r>0\), this is a product of the structure of the equation. The fixed points are \(x=0\) and \(x=\pm\sqrt{r}\), given by solving \(f^\prime(x)=0\) and the derivative being \(f^\prime(x)=r-3x^2\), and these points results in \(f^\prime(0)=r\) and \(f^\prime(\pm\sqrt{r})=-2r\). Given that a system is stable if \(f^\prime(x)<0\), more about this in section 6.3 Stability properties of an equilibrium, \(f^\prime(r)\) is stable for \(r<0\) and unstable for \(r>0\). As such, when considering \(r>0\), the point \(x=0\) is unstable and the two new stable points emerges at \(x=\pm\sqrt{r}\).

Hopf bifurcation

This bifurcation type is an example of an two-dimensional bifurcation and like the pitchfork bifurcation there is a supercritical and subcritical subtype. In the case of the supercritical type, an unstable equilibrium point with and a stable limit cycle emerges from a stable equilibrium when crossing the bifurcation point. For the subcritical type, a stable equilibrium with an unstable limit cycle transforms to an unstable equilibrium when crossing the bifurcation point. To illustrate, we consider the polar coordinates:

\[ \begin{bmatrix} \frac{dr}{dt} = (\mu - r^2)*r \\ \frac{d\theta}{dt} = w \end{bmatrix} \]

Here, \(r(t)\) describes the amplitude of an oscillation and \(\theta(t)\) describes the angular position. Skipping the analytical analysis (see task 6.3), this system behaves accordingly:

  • \(\mu>0\), \(r(t)\) has an unstable equilibrium point (at \(r=0\)) and a stable equilibrium point (at \(r=\sqrt{\mu}\)) with a stable limit cycle with radius \(\sqrt{\mu}\).
  • \(\mu<0\), \(r(t)\) has only a stable equilibrium point (at \(r=0\)).

This is then an example of a supercritical hopf bifurcation, if we instead change the sign of \(\frac{dr}{dt}\)

\[ \begin{bmatrix} \frac{dr}{dt} = -(\mu - r^2)*r \\ \frac{d\theta}{dt} = w \end{bmatrix} \]

We get a subcritical type as we reverse the stability of the equilibrium points of \(r(t)\):

  • \(\mu>0\), \(r(t)\) has a stable equilibrium point (at \(r=0\)) with an unstable limit cycle with radius \(\sqrt{\mu}\) and an unstable equilibrium point (at \(r=\sqrt{\mu}\)).

For clarification, a stable limit cycle indicates that all trajectories in the system is drawn to the limit cycle as \(t\to\infty\). Whereas as trajectories is instead drawn to an unstable limit cycle when \(t\to-\infty\). We can thus illustrate these behaviours with the diagram below, where the trajectory is given in red, a stable limit cycle in solid blue, and an unstable limit cycle in dashed blue:

Hopf bifurcation

These are just a few examples of bifurcations in mathematical models. Bifurcation analysis provides a powerful tool for understanding the behaviour and dynamics of complex biological systems, helping researchers gain insights into the underlying mechanisms and predict system responses to changes in parameters or external stimuli.

Since bifurcations occur when the equations governing the system undergo a qualitative change in their solutions, due to changes in the model's parameters, modellers must foresee and identify these bifurcations. To do this one can perform a bifurcation analysis on the chosen system and this is what we will do in this section. Before going into this analysis, you can read a bit more about it below and familiarize yourself with the model we will analyze, the Hodgkin-Huxley model.

Bifurcation analysis

Bifurcation analysis is a method used to study the qualitative changes in the behaviour of a system as a parameter or set of parameters is varied. It provides insights into the different states and dynamics that a system can exhibit, such as the presence of equilibrium points, periodic orbits, or chaotic behaviour.

The goal of bifurcation analysis is to identify critical parameter values at which qualitative changes occur and to understand how these changes arise. By examining the bifurcations, one can determine the stability of equilibrium points, the existence of periodic solutions, and the emergence of complex behaviour.

Here are some steps typically involved in bifurcation analysis:

  • Define the mathematical model: Start with a set of differential equations or discrete-time equations that describe the system under study. These equations should capture the dynamics and interactions of the variables in the system.

  • Identify relevant parameters: Determine the parameters in the model that can be varied to investigate the system's behaviour. These parameters could represent physical constants, external inputs, or other factors that influence the system's dynamics.

  • Analyze equilibria: Determine the equilibrium points of the system by setting the derivative equations to zero. This involves solving a set of algebraic equations. For each equilibrium point, calculate its stability properties by analysing the eigenvalues of the linearized system.

  • Perform parameter continuation: Choose a parameter of interest and vary its value over a specified range. By using numerical techniques, such as the pseudo-arclength continuation method or Newton's method, you can track the equilibrium points and their stability as the parameter is varied.

  • Detect bifurcation points: Monitor the behaviour of the equilibrium points as the parameter is varied. Look for critical parameter values where qualitative changes occur, such as the appearance or disappearance of equilibrium points or the transition from stability to instability. These critical points are bifurcation points.

  • Classify bifurcations: Determine the type of bifurcation occurring at each critical point. This involves analysing the eigenvalues, normal forms, and other characteristics of the system near the bifurcation point. Common bifurcation types include saddle-node, pitchfork, Hopf, and Turing bifurcations, among others.

  • Explore parameter space: Systematically explore the parameter space to identify regions with different dynamic behaviours. Plot bifurcation diagrams, which display the stable and unstable equilibrium points or periodic orbits as a function of the parameter(s). These diagrams provide an overview of the system's behaviour and the regions where different dynamics occur.

Bifurcation analysis is a powerful tool for understanding the behaviour of complex systems and how their dynamics change with varying parameters. It allows modellers to uncover critical parameter values and gain insights into the underlying mechanisms that give rise to different qualitative behaviours.

Hodgkin-Huxley model

The Hodgkin-Huxley model is a mathematical model that describes the electrical activity of neurons and is widely used to understand the generation and propagation of action potentials.

The Hodgkin-Huxley model represents the neuronal membrane as a circuit of ion channels that allow the flow of specific ions, such as sodium (Na+) and potassium (K+), across the cell membrane. The changes in ion flow across the membrane lead to changes in the membrane potential and the generation of action potentials. The model consists of a set of coupled differential equations that describe the behaviour of the ion channels and their interactions. The main variables in the model are:

  • V: Membrane potential of the neuron.
  • m, h, and n: Activation variables that represent the opening and closing of specific ion channels (sodium and potassium) in response to changes in voltage.

The model's equations can be simplified as follows:

dV/dt = (1/C) * (I - gNa * m^3 * h * (V - ENa) - gK * n^4 * (V - EK) - gL * (V - EL))
dm/dt = αm * (1 - m) - βm * m
dh/dt = αh * (1 - h) - βh * h
dn/dt = αn * (1 - n) - βn * n

Here, \(C\) is the membrane capacitance, \(I\) represents the applied current, \(gNa\) and \(gK\) are the conductances of sodium and potassium channels, respectively, \(ENa\) and \(EK\) are the reversal potentials for sodium and potassium ions, \(gL\) and \(EL\) represent the conductance and reversal potential of the leakage channels, respectively. The terms \(α\) and \(β\) in the equations represent voltage-dependent rate constants that control the opening and closing of ion channels. These rate constants are functions of the membrane potential (\(V\)) and determine the kinetics of the ion channels, we will not go into the numerical expressions that formulate these rate constants.

The Hodgkin-Huxley model describes how changes in membrane potential and the interactions between ion channels result in the generation and propagation of action potentials. When the membrane potential reaches a certain threshold, the sodium channels open, allowing an influx of sodium ions and causing a rapid depolarization of the membrane. This depolarization triggers a positive feedback loop, leading to the rapid rise of the membrane potential, known as the upstroke of the action potential.

After the upstroke, the potassium channels open, allowing an efflux of potassium ions, which repolarizes the membrane and restores it to its resting state. The closing of sodium channels and the opening/closing of potassium channels are controlled by the activation variables (\(m\), \(h\), and \(n\)), which are influenced by voltage-dependent rate constants. The model can be reformulated be a bit easier to understand, se below:

dV/dt = (1/C) * (I + INa + IK - IL) 
dm/dt = αm * (1 - m) - βm * m
dh/dt = αh * (1 - h) - βh * h
dn/dt = αn * (1 - n) - βn * n

where,
INa = - gNa * m^3 * h * (V - ENa)
IK = - gK * n^4 * (V - EK)
IL = gL * (V - EL)

In this format, it is easy to understand that the currents from the sodium (\(INa\)) and potassium (\(IK\)) channels increases the membrane potential (\(V\)) and the leakage current (\(IL\)) decreases the membrane potential. Where, the parameters \(gNa\), \(gK\), and \(gL\) are scaling these currents. You will focus on the currents (\(INa\), \(IK\), \(IL\)) and the potential (\(V\)) in this lab and the overall interaction between the membrane potential and the three currents. It is therefore enough to remember that the activation variables \(m\), \(n\), and \(h\) describe the activation (\(m\) and \(n\)) of the sodium and potassium channels respectively, while (\(h\)) describe the deactivation of the sodium channels.

The Hodgkin-Huxley model accurately reproduces the behaviour of real neurons and has been instrumental in understanding the mechanisms underlying the generation and propagation of action potentials. It provides valuable insights into various neurological phenomena and has paved the way for further advancements in neuroscience.

Now then, we will follow the introduced steps for the bifurcation analysis and identify the model behaviour over some parts of the parameter space.

6. Bifurcation analysis of the Hodgkin-Huxley model

The first step has already been decided for you, we will analyze the Hodgkin-Huxley model. You have been provided with a prepared script BifurcationAnalysis to assist you throughout the analysis and you can uncomment code sections as we advance throughout the steps.

The first task for you is to observe the model structure and try to get an idea of how the parameters influence the model behaviour.

6-2. Identify relevant parameters

Observing the model structure, we can identify that the variable parameters are

  • the conductances of the sodium and potassium channels, \(gNa\) and \(gK\)
  • the reversal potentials, \(ENa\) and \(EK\)
  • the leak conductance, \(gL\), with its corresponding reversal potential \(EL\)

To be able to know that the changes in model behaviour come from changes of the parameter values we must first find an equilibrium. Remember that the currents from the sodium (\(INa\)) and potassium (\(IK\)) channels increases the membrane potential (\(V\)) and the leakage current (\(IL\)) decreases the membrane potential. Meaning, that INa+IK-IL=0 must be true for the model to be in equilibrium (if no external stimulus is applied, i.e. I=0).

6-3. Analyze equilibriums

As it can be a bit tricky to find all equilibriums, the example script BifurcationAnalysis is prepared to do so for you. The script will also plot the equilibrium for the model, relating them to the different currents. This section does quite a few things, to shortly summarize them:

  1. First, the equilibriums are calculated. Scipy's fsolve is utilized to algebraically solve the equilibriums for the Hodgkin Huxley model. As we learnt previously in the Domains of attraction section, it can be difficult to find the equilibriums. In this case, there are 3 equilibriums and the function solveEQ gives fsolve three different initial guesses to find all three of them for you.
  2. The membrane potential is simulated and plotted for all three equilibriums to verify that the model is in steady state (d/dt = 0).

  3. The stability of the equilibriums are investigated. Using Scipy's eig function, the eigenvalues can be calculated from the Jacobian (a 4x4 matrix with each model state derived against themselves). The math has been done for you and the Jacobian is already defined. You can find resourced of this process below.

    Stability properties of an equilibrium

    To numerically analyze the stability of an equilibrium we would first need to linearize the system, we can fo this with the help of the jacobian matrix, J.

    \[ J= \begin{bmatrix} J_{11}& J_{12}& J_{13}& J_{14}\\ J_{21}& J_{22}& J_{23}& J_{24}\\ J_{31}& J_{32}& J_{33}& J_{34}\\ J_{41}& J_{42}& J_{43}& J_{44} \end{bmatrix} = \begin{bmatrix} \frac{df_V}{dV}& \frac{df_V}{dm}& \frac{df_V}{dn}& \frac{df_V}{dh}\\ \frac{df_m}{dV}& \frac{df_m}{dm}& \frac{df_m}{dn}& \frac{df_m}{dh}\\ \frac{df_n}{dV}& \frac{df_n}{dm}& \frac{df_n}{dn}& \frac{df_n}{dh}\\ \frac{df_h}{dV}& \frac{df_h}{dm}& \frac{df_h}{dn}& \frac{df_h}{dh}\\ \end{bmatrix} \]

    In the case of the states \(m\), \(n\), and \(h\) they are not dependent on each other, meaning that: \(\frac{df_m}{dn}\), \(\frac{df_m}{dh}\), \(\frac{df_n}{dm}\), \(\frac{df_n}{dh}\), \(\frac{df_h}{dm}\), and \(\frac{df_h}{dn}\) equals 0. The cellmembrane potential state, \(V\), are straight forward to derivate and the exponent are just reduced for \(\frac{df_V}{dm}\), \(\frac{df_V}{dn}\), and \(\frac{df_V}{dh}\). For \(\frac{df_V}{dV}\), the constant part (\(V-E_i\)) of the currents are just reduced to \(*1\). The derivatives of \(\frac{df_m}{dm}\), \(\frac{df_n}{dn}\), and \(\frac{df_h}{dh}\) are more complex to arrive at and you need to make use of chain rule and at some places substitution. The final Jacobian looks like:

    \[ J= \begin{bmatrix} \frac{df_V}{dV}& \frac{df_V}{dm}& \frac{df_V}{dn}& \frac{df_V}{dh}\\ \frac{df_m}{dV}& \frac{df_m}{dm}& 0& 0 \\ \frac{df_n}{dV}& 0& \frac{df_n}{dn}& 0 \\ \frac{df_h}{dV}& 0& 0& \frac{df_h}{dh} \end{bmatrix} \]

    To find the eigenvalues, we need to solve the characteristic equation \(det(J−\lambda I)=0\). Here, \(I\) is the identity matrix and \(\lambda\) represents the eigenvalues.

    \(J−\lambda I\) is formulated as:

    \[ J-\lambda I = \begin{bmatrix} \frac{df_V}{dV}−\lambda& \frac{df_V}{dm}& \frac{df_V}{dn}& \frac{df_V}{dh} \\ \frac{df_m}{dV} & \frac{df_m}{dm}−\lambda& 0& 0 \\ \frac{df_n}{dV} & 0& \frac{df_n}{dn}−\lambda& 0 \\ \frac{df_h}{dV} & 0& 0& \frac{df_h}{dh}−\lambda \end{bmatrix} \]

    Solving the determinant for \(det(J-\lambda I)=0\) gives us four eigenvalues, by analysing these one can interpret how stable a system is and how it will behave. For instance, some of the things to look out for are:

    • The real part of the eigenvalues are all negative (stable phase) - indicating that the system is stable.
    • Any real part of the eigenvalues are positive (unstable phase) - indicating that the system is unstable.
    • The imaginary part are not zero - indicating the system is subject to oscillations.
    • The eigenvalue are equal to zero (center space) - indicating that this equilibrium point is a non-hyperbolic point (often bifurcation point).

    For our equilibrium points we can observe that:

    • Equilibrium A is stable (all eigenvalues is negative).
    • Equilibrium B is unstable (as one eigenvalue is positive).
    • Equilibrium C is unstable (as multiple eigenvalues are positive) and is subject to oscillations (imaginary part is not zero).

    To further get an understanding of the behaviours, we will continue to analyse the equilibriums graphically.

  4. We simulate the model and record the currents over a fixed range of voltages, -85 to -40 mV, which is of physiological interest. Here, we get the behaviour of the currents over the voltage range.

    cvHandleFailure

    As we are simulating over the currents, -85 to -40 mV, we will cross the bifurcation point and as such a simulation will crash and produce the error [ERROR][rank 0][sundials-7.0.0\src\cvode\cvode.c:3698][cvHandleFailure] At t = 0 and h = 1.3093e-10, the corrector convergence test failed repeatedly or with |h| = hmin.. This is expected and do not indicate that anything is wrong.

  5. The script plots the sodium current, potassium current, combined sodium and potassium current, leakage current, and the three equilibrium points into one figure EQ currents to illustrate their behaviours as the membrane potential varies.

In the acquired figure EQ currents, we can see that the leakage current is intersected three times, at the locations for the equilibriums. Equilibrium A refers to the resting state and is a stable equilibrium. Equilibrium C refers to the maximum action potential value the model can achieve during the process of depolarization. Equilibrium B is a transient state between the resting state and the depolarization state. As we change the parameters of the model, equilibrium B will shift but remain stable. The model state at equilibrium C, on the other hand, is not approaching stability. As we offset the model, by changing the parameters, the model state at C can cause a change of topological structure and behaviour. Moving forward, we will further analyse the third equilibrium (C).

With the knowledge of which points in the parameter space that the model is stable, we can now proceed with the parameter continuation analysis. Our goal here is to learn how the model behaves over a physiological range of parameter values and potentially a broader range, since we would use the model within these ranges.

6-4. Perform parameter continuation

Now let us investigate how the leakage conductance parameter, gL, affects the model. The script is prepared to iterate over a range of \(gL\) values and then to plot the basal, steady state, membrane potential against the value of \(gL\). Try different ranges of \(gL\), however stay within 0-1.0 mS/cm^2.

Observing the behaviour of the basal membrane potential over the range of \(gL\) values. Is there a distinct change of behaviour for a \(gL\) value?

6-5. Detect bifurcation points

With support from the figure you got in the previous section, select a few values of \(gL\) and add them to gL_values under Step 5. The script will plot the cell membrane behaviour over time, with the initial guess being equilibrium C. We can then see how the cell membrane potential tunes into the resting potential for the \(gL\) values you entered in the gL_values. For which value of \(gL\) do you find a change of this behaviour? Note that: - You might need to change the resolution (the step size) of the \(gL\) values. - As we are deriving the equilibrium point numerically, we will not be able to find the exact value of \(gL\) that causes the bifurcation point. A resolution of 2 decimals (.xx) is enough. - The oscillations might not always have had enough simulation time to transition into the equilibrium value. You could increase the simulation time if you want to or you focus on the general behaviour/trend of the oscillations.

Equilibrium C - Unstable

Remember that equilibrium C was unstable, which allows the model to find a resting membrane potential that is not found for any other value of \(gL\) in the same region. This is not a bifurcation point.

When a deviant behaviour is found we want to identify what happens. Different types of bifurcations have distinct behaviours and can be identified from the change in model behaviour.

6-6. Classify bifurcations

Going back to definitions of different bifurcations given in the introduction, three bifurcations were introduced.

  • Saddle-node bifurcation: In this bifurcation, an equilibrium point collides with and disappears from the system as a parameter is varied. It leads to the creation or destruction of stable equilibrium points.

  • Pitchfork bifurcation: In this bifurcation, an equilibrium point undergoes a symmetry-breaking transition, resulting in the emergence of new equilibrium points (supercritical) or the disappearance of equilibrium points (subcritical).

  • Hopf bifurcation: In this bifurcation, a stable equilibrium point loses stability, and a stable limit cycle emerges as a parameter is varied. It leads to the appearance of periodic oscillations in the system as \(t\to\infty\) (supercritical). If instead the oscillations occur when \(t\to-\infty\) it is the subcritical variant.

What have you found from observing the behaviour in the figures from part 3, part 4, and part 5? Note that the figures share common axises, although the membrane potential is sometimes on the y-axis and sometimes on the x-axis. Can the behaviour of the bifurcation point be described by one of these bifurcations?

To summarize the observed behaviour in the previous sections, we stated in the unstable equilibrium C. By increasing the value of \(gL\) the leakage current increases, which causes the membrane potential to oscillate increasingly aggressively before returning to a resting state. Eventually, at the bifurcation point, the model will not stop oscillating and never return to a resting potential. As we further increase the leakage conductance, the model leaves equilibrium C and falls back to equilibrium A, at a resting potential approaching -81 mv. In this case, we have an unstable equilibrium point, C, with a stable limit cycle (the oscillating behaviour) surrounding it - making this bifurcation point a Hopf bifurcation, more specifically a Hopf bifurcation of the type supercritical.

The last part of the analysis is to explore the parameter space and familiarize us with the total model behaviour.

6-7. Explore the parameter space

Part 7 in the script is set up to plot a 3D plot describing the transition of the states controlling the opening of the ion channels, m - the opening state of Na+, h - the closing state of Na+, and n - the opening state for K+. Plot the behaviour for values of \(gL\) i) lower than the bifurcation value, ii) the bifurcation value, and iii) over the bifurcation value (ignore the unstable point for equilibrium C). You should find three different behaviours for these three cases. Note that the red dot marks the starting point of the transition.

  1. for \(gL\) values under the bifurcation point the model quickly find a stable position and the states will stabilize into a point - indicating a stable position.
  2. for \(gL\) values equal to the bifurcation point the model will showcase a transition into stable and constant oscillating behaviour
  3. for \(gL\) values above the bifurcation point to model quickly find a stable position, different from the one in case 1, and the states will stabilize into a point - indicating a stable position.

As the model transition towards its behaviour, the circulating motions can be explained by the oscillating behaviour, before settling at a resting position, of the cell membrane potential you observed in part 5. Additionally, at the bifurcation point the oscillating behaviour never seized, which is illustrated here by the repeating loops.

Now we have performed a thorough analysis for the parameter governing leakage conductance. To use this model, and other nonlinear model, with confidence - you need to be aware of its capabilities and faulty behaviours. The unstable equilibrium C represent a model state where the neuron would not leave the hyperpolarisation peak, something that never happens. Neither is the almost immediate transition over the bifurcation point a good physiological description. These behaviours are introduced by the model design. This exercise has hopefully given you some insights regarding bifurcation analysis, if you want to look at some more thorough work you can look at the following articles i) Analysis and control of the bifurcation of Hodgkin–Huxley model, ii) Effects of Maximal Sodium and Potassium Conductance on the Stability of Hodgkin-Huxley Model.

Now, perform the analysis again for the parameter governing the conductance of the sodium channel, gNa. You can easily repurpose the BifurcationAnalysis script for this, ideally you make a different script for this.

7. Bifurcation analysis of sodium channel

Repeat the analysis steps that you did for the leakage conductance parameter but perform the analysis on the sodium conductance parameter instead. You can start at part 4, since the equilibriums will be the same, and explore the behaviour from equilibrium C. Is the model also sensitive to changes of \(gNa\)? Is there a bifurcation point for \(gNa\)? A reasonable range to explore \(gNa\) over is 0-500 mS/cm^2. Remember that you cannot find the exact value of \(gNa\) using a graphical approach, a resolution of 1 decimal (.x) is enough.

What value of gNa causes a bifurcation?

What value of gNa causes a bifurcation?

Questions for lab 1B

In your lab report, you should answer the following questions:
  • Describe what the "transient time" is.
  • What are the benefits of a stable equilibrium?
  • How can non-stable equilibrium influence our models?
  • How are semi-stable equilibriums defined?
  • Describe a type of bifurcation (of your own choice).
  • What is the general strategy of a bifurcation analysis?
  • Is our example of Hodgkin-Huxley sensitive to changes of the sodium conductance parameter? Are there bifurcation points? Please motivate your answer with images of your results.