TBMT42 Lab 1.2: 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 behaviors vs stationary model behaviors and how they can change over time. Secondly, we will investigate the different types of equilibriums there are and how domains of attraction influences 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 behaviors of models and is very useful when working with nonlinear systems, which might have unforeseen (or hard to predict) behaviors 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 Matlab. You can download the example files here. You are free to use any software, however, if you don't use Matlab you will need to redo the pre-made scripts yourself. For installation and setup of, see Matlab installation and Matlab packages.

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 behaviors, 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 behavior. Nonlinear systems can exhibit phenomena such as oscillations, chaos, bifurcations, and stability transitions, which are critical for understanding the behavior 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 behavior 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 behavior of neural networks.

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

Transient stationary behavior¶

Concept of transient and stationary behavior

Firstly, let's delve into the concept of transient and stationary behavior. 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 behavior. 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 behaviors.

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 behavior

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 TransientStationaryBehavior.m 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 behavior

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 state in the model. This is one way of implementing time dependent effects, where a constant state (derivate=0) can be used to set the value of the input effect from the simulation call. This allows us to simulate the 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 optimizing call. Below you see the state implementation.

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

where meal_influx = 0.1 if 0 < t <= t_on or meal_influx = 0 otherwise.

This has been implemented in the TransientStationaryBehavior.m as a simulation switch, where the model is simulated with meal_influx = 0.1 for timeDuringMeal minutes and then additionally simulated with meal_influx = 0 for timeAfterMeal minutes. The switch is implemented in the Utilities/SimulateWithInput function, which takes the model’s name, parameter values, the first time period, state ID, the first input value, second time period, and second input value as inputs. The function then returns a concatenated structure with the simulation results from both simulations.

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 meal_end variable in the TransientStationaryBehavior 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 behavior, 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

Return to the timeSimulationEnd variable and change the post-meal duration for the simulation (the value must be greater than the meal_end value), also remember to uncomment the Transient behavior 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 (theta1 and theta2). 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 behavior.
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 behavior 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 behavior 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 modeling of the COVID-19 and its dynamics.

The equations for our SEIR model are as follows:

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

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 behavior 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 is some question regarding this behavior: - 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 observe the equilibrium point for the logistic growth model. This model has two equilibrium points, which you can read more about below. In essence, one equilibrium is unstable and can only be found from one starting position. You can read more about the model below.

The logistic growth model

We start by looking at the simple logistic growth model, given by

d/dt(P) = r * P * (1 - P/K)

Here, $P$ represents the population size, $t$ represents time, $r$ is the intrinsic growth rate, and $K4$ is the carrying capacity of the environment.

5. Find the unstable equilibrium

The function UnstableLogisticGrowthModel is prepared for you to illustrate the equilibrium points and how the initial population size P(0) affects this. The script will give you a figure Logistic growth model equilibriums that plot the equilibrium values for the population P for three different initial guesses. Analyze the figure and find the initial guess that gives the unstable equilibrium.

What starting value finds the unstable equilibrium

For what starting value P(0) is an unique, unstable, equilibrium found?

Explanation of the behavior

The unstable equilibrium point represents either extinction or the absence of a population. The other equilibrium point is found at P = K, where K is the carrying capacity of the environment, i.e. where the population size stabilizes and growth stops.

If the population size is perturbed slightly away from the equilibrium point ´P = K´, either above or below the carrying capacity, the dynamics of the system will return to this equilibrium. While the unstable equilibrium point for this model only can be found if there is no existing population to being with. As soon as some individuals are introduced to the population, the population will move towards the carrying capacity. In practice this feature of the model is not that useful, since we don't need a model to say that a population cannot grow if there are no individuals.

Moving on, we will now do some more work with the FitzHugh-Nagumo model, a more complex nonlinear 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:

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.

6. 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 behavior 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 behavior if that behavior has a constant period time. What happens with the model behavior as we approach the unstable equilibrium?

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 behavior 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 behavior 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 behavior exhibited by the system.

Bifurcations are important in biological modeling because they can help explain the observed transitions and changes in behavior 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.
Pitchfork bifurcation: In this bifurcation, an equilibrium point undergoes a symmetry-breaking transition, resulting in the emergence of two new equilibrium points. It can lead to the creation or destruction of stable or unstable equilibrium points.
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.
Turing bifurcation: This bifurcation occurs in spatially extended systems and is associated with the formation of spatial patterns. It can lead to the emergence of complex spatial structures and self-organization in biological systems.

These are just a few examples of bifurcations in mathematical models. Bifurcation analysis provides a powerful tool for understanding the behavior 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, modelers 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 behavior 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 behavior.

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 behavior.

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 behavior. 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 analyzing 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 behavior 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 analyzing 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 behaviors. 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 behavior and the regions where different dynamics occur.

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

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 behavior 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 $gNa4, \(gK$, and 4gL\) 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 $h4 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 behavior 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 behavior over some parts of the parameter space.

7. 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 behavior.

7-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, 4gL$, with its corresponding reversal potential $EL$

To be able to know that the changes in model behavior 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).

7-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:

First, the equilibriums are calculated. Matlab's vpasolve 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 here and solveEQHH gives vpasolve three different initial guesses to find all three of them for you.
Using an equilibrium solution we solve the currents for a fixed range of voltage, -85 to -40 V, which is of physiological interest. Here, we get the behavior of the currents over the voltage range.
We also simulate the model for every equilibrium to get the currents in the equilibrium points.
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 behaviors as the membrane potential varies.
Lastly, an additional figure Steady state of EQ A, B, and C is plotted to showcase that the model is in steady state for these equilibriums and that $\frac{dV}{dt}=0$ holds true.

The behavior of the currents

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 behavior.

These equilibriums could be solved numerically, which also would allow us to analyze the behavior numerically. We will not pursuit this path, instead we will focus on studying the behavior graphically. If you are interested in how one would go about doing this numerically, a summarized description is available 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}{dh}& \frac{df_V}{dn}\\ \frac{df_m}{dV}& \frac{df_m}{dm}& \frac{df_m}{dh}& \frac{df_m}{dn}\\ \frac{df_h}{dV}& \frac{df_h}{dm}& \frac{df_h}{dh}& \frac{df_h}{dn}\\ \frac{df_n}{dV}& \frac{df_n}{dm}& \frac{df_n}{dh}& \frac{df_n}{dn} \end{bmatrix} = \begin{bmatrix} \frac{df_V}{dV}& \frac{df_V}{dm}& \frac{df_V}{dh}& \frac{df_V}{dn}\\ \frac{df_m}{dV}& \frac{df_m}{dm}& 0& 0\\ \frac{df_h}{dV}& 0& \frac{df_h}{dh}& 0\\ \frac{df_n}{dV}& 0& 0& \frac{df_n}{dn} \end{bmatrix} \]

We can formulate our system linearized with the in the following way.

\[ \begin{gather*} \frac{dV}{dt} = J_{11}*V + J_{12}*m + J_{13}*h + J_{14}*n \\ \frac{dm}{dt} = J_{21}*V + J_{22}*m \\ \frac{dh}{dt} = J_{31}*V + J_{33}*h \\ \frac{dn}{dt} = J_{41}*V + J_{44}*n \end{gather*} \]

Here, the characteristic equation is formulated as:

\[ \lambda^4 + a\lambda^3 + b\lambda^2 + c\lambda + d = 0 \\ \]

where,

\[ \begin{gather*} a = -(J11 + J22 + J33 + J44) \\ b = J11*(J22 + J33 + J44) + J22*(J33 + J44) + J33*J44 - J12*J21 - J13*J31 - J14*J41 \\ c = J12*J21*(J33 + J44) + J13*J31*(J22 + J44) + J14*J41*(J22 + J33) - J11*J22*(J33 + J44) - (J11 + J22)*J33*J44 \\ d = J11*J22*J33*J44 - J12*J21*J33*J44 - J13*J22*J31*J44 - J14*J22*J33*J41 \\ \end{gather*} \]

According to the Routh-Hurwitz criterion the roots of the linearized system are all negative if.

\[ \begin{gather*} a>0\\ a*b>c\\ a*b*c>c^2 + a^2*d\\ d>0\\ \end{gather*} \]

Meaning that the system is stable. This can further be observed by solving the roots from of the system $(J-\lambda I)$, where

\[ \lambda I = \begin{bmatrix} \lambda I& 0& 0& 0\\ 0& \lambda I& 0& 0\\ 0& 0& \lambda I& 0\\ 0& 0& 0& \lambda I \end{bmatrix} \]

Giving us

\[ (J- \lambda I) = \begin{bmatrix} J11-\lambda I& J12& J13& J14\\ J21& J22-\lambda I& 0& 0\\ J31& 0& J33-\lambda I& 0\\ J41& 0& 0& J44-\lambda I \end{bmatrix} \]

Solving the determinant for $det(J-\lambda I)=0$ gives us four eigenvalues, by analyzing these one could determine the bifurcation points. For instance, some of the things to look out for to determine a bifurcation point are: - The eigenvalues change sign. - The eigenvalues cross the imaginary axis. - The eigenvalues becomes zero.

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.

7-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 behavior of the basal membrane potential over the range of $gL$ values. Is there a distinct change of behavior for a $gL$ value?

7-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 behavior 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 behavior? 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 behavior/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 behavior is found we want to identify what happens. Different types of bifurcations have distinct behaviors and can be identified from the change in model behavior.

7-6. Classify bifurcations

Going back to definitions of different bifurcations given in the introduction, four 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 two new equilibrium points. It can lead to the creation or destruction of stable or unstable equilibrium points.
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.
Turing bifurcation: This bifurcation occurs in spatially extended systems and is associated with the formation of spatial patterns. It can lead to the emergence of complex spatial structures and self-organization in biological systems.

What how you found from observing the behavior in the figures from part 4 and part 5? Can the behavior of the bifurcation point be described by one of these bifurcations?

To summarize the observed behavior 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 behavior) 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 behavior.

7-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 behavior 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 behaviors for these three cases.

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.
for $gL$ values equal to the bifurcation point the model will showcase a transition into stable and constant oscillating behavior
for $gL$ values above the bifurcation point to model quickly find a stable position, different from the the one in case 1, and the states will stabilize into a point - indicating a stable position.

As the model transition towards its behavior, the circulating motions can be explained by the oscillating behavior, before settling at a resting position, of the cell membrane potential you observed in part 5. Additionally, at the bifurcation point the oscillating behavior 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 behaviors. The unstable equilibrium C represent a model state where the neuron would not leave the hyperpolarization peak, something that never happens. Neither is the almost immediate transition over the bifurcation point a good physiological description. These behaviors 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.

8. 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 behavior 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?

Questions for lab 1.2¶

In your lab report, you should answer the following questions:

What is the transient time?
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 potassium conductance parameter? Are there bifurcation points?