It is never easy or cheap to carry out biological research, without a laboratory, a microscope, etc. In general, practical constraints makes challenging to design practical and portable demonstrations and outreach activities to illustrate key findings and concepts. On the other hand, using mathematical models based on physical and chemical principles allow us to formulate the theoretical frameworks which can be investigated using the machinery of dynamical systems and computer simulations, which lead to hypotheses testing, and/or predictions about the mechanisms under investigation, again, unless the target audience for a demonstration or talk is not versed with such concepts, it can be tricky to make translate any results based only on equations and static plots.
In this post, I present a set of examples which, couple experiments, maths and programming nicely in a way which can be presented to biologists, physicists and general audiences to illustrate how using in the theory, experiments alongside with computers can be used to illustrate ideas is, genetics, cardiac dynamics, neurosciences alike using a simple mechanical system.
Again, for this project all the simulations I show here, had all the numerical codes written in C and compiled using emscripten to be used in the simulations which are all built using my favorite JS library three.js.
Excitable Systems - Shishi-odoshi.#
I would like to put under consideration the oscillator illustrated below:
Prettier versions of this device can be found in some Japanese gardens and are called “deer-frigthening” or “boar-frigthening” and according to Wikipedia that has been the traditional use of these. Beyond beauty and elegance, these devices are interesting because these show a set of features which are shared by some very important biological oscillators:
The system operates due to an external stimulus.
There exists an activation threshold, which needs to be reached before the system activates.
There is a refractory period. Which is a time interval during which the system remains inactive even if it is under stimulation.
The little movie shown below (which I borrowed from Flavio Fenton’s old website) shows a myocyte being stimulated at a single location. The stimulus location releases calcium which then diffuses through the cell, activating more release locations sites and forming a wave. The stimulus is continuous, yet the behaviour of the myocyte is that of a beating cell.
Propagating calcium wave producing the contraction of a ventricular cell. Copyright: Lee & Fenton
Models and Arrhythmias.#
(A): Human atrial cell compartment model schematics (B): The RyRs states Markov model with the four configurations O: open, C: closed, and Inactive states I1 and I2.
The diagram above shows: (A) the way all cell compartments are represented and connected through pumps and gates, connecting several inner and outer domains of the cell. (B) The complex responsible for the calcium cycling in these type of cells are the Ryanodine receptors (RyR2) which in this model can be found in four configurations or states \( (I_1, I_2, O,C) \). The reciprocal of the transition rate \( \tau_r = \frac{1}{k_B} \) corresponds to the recovery from inactivation time.
Left: Simplified cell compartment model schematics. Right: The RyRs states Markov model with the four configurations O: open, C: closed, and Inactive states I1 and I2.
The translation of the dynamics sketched in the diagrams into mathematical models can be as complex and detailed as we require (see [1]). In this post however I would like to use the model presented in [2], which serves very well to illustrate what we want in this post.
$$\dot{c}_{j} = I_{s}(t)+g_{rel}c_{rel}\frac{k_ac_j^2}{k_{om}+k_ac_j^2}(c_{SR}-c_j) - (c_j-c_o)/\tau_{diff}\quad(1)$$
$$\dot{c}_{rel} = k_{im}(1-c_{rel}) - k_ic_jc_{rel}\quad(2)$$
(1) and (2) constitute a minimal model for calcium release incorporating the release dynamics. Notice that the first term in (1)’s rhs is the external stimulus \(I_{CaL}(t)=I_{s}(t)\), Which is provided by the pacemaker cells in the heart or externally in experiments. These cells provide a periodic stimulus \(T_s\). This is where things get interesting.
A healthy cell will normally release Calcium periodically, driven by the stimulus if the recovery time of the ion release channels is slower compared to the stimulus. An example is illustrated below. Where \(T_s=400\ ms\) and \(T_r = 200\ ms \)
RyR2 malfunctions, Alternans and the origin of Fibrillation.#
If on the other hand the release is slower due to a malfunction of the calcium release units, which might lead to a slow recovery, the periodic response alternates between a high/low beat to beat calcium release regime. These calcium alternans are a known mechanism for inducing cardiac arrhythmias such as atrial fibrillation. Below I show results for the case of \( T_s=400\ ms \) and \(T_r = 320\ ms \) .
Stimulus-Recovery relationship and Interactive Solver.#
Bifurcation plot: Maximum amplitudes of Cj vs τrec It can be seen which intervals present regular oscillations, 2:1 alternans as well as 3:1 alternans.
Back to the seesaw#
The complex alternating behaviour of the cardiac cells occurs if the cell possess dysfunctional calcium release units due to a damaged cell, or an underlying genetic condition. Mathematically this is encoded in the value of \( k_{im} \).
To investigate the origin of these dynamics, it is easy to build a japanese seesaw with bits and bobs of Lego and K’Nex sets (The best ever construction set available? Certainly the most versatile!). And make the following observations about the system parameters.

Alongside these condition we need to incorporate the stimulus, in this case by adding mass to the mass attached in the right end of the rod. A very easy way to do this is to use the conditions:
$$m_2=m_{2o} + \alpha t$$
If the angle \(\Phi \geq \Phi^{\dagger} \). Where \(\Phi^{{\dagger}}\) is a parameter stating the range of action if the stimulus. And the condition \( m_2(t') = \gamma m_2(t) \) if \( \Phi = \Phi_2\). Where \(\gamma \in [0,1] \) is a parameter stating the fraction of water released if we reached the boundary condition at \( \Phi_2 \). This parameter is the analogous to the parameter \( k_{im}\), As it regulates the recovery time.
We can see that there are a few parameters which are important, for instance the threshold mass \(m_1 \) which needs to be surpassed by \(m_2\) to move the centre of mass from the quadrant satisfying the constraint \(\Phi_1\) into motion. Another two important parameters are the current \(\alpha\) as well as the fraction of water released mentioned above. We will return to this at the end of the post. For now, let us see the system in action!
Experiments.#
I

The Japanese seesaw toy: (Left and Middle left) The device. (Middle Right and Right) The replaceable release mechanism: by using different configurations we can modify the release of fluid, therefore the refractory period.
The seesaw in action! This is the periodic behaviour we can expect, if the release in almost complete and the refiill is not very fast, so the system fully recovers its initial position.
The seesaw in action! This is the WT expected periodic behaviour if the release in almost total and the refill is not very fast, allowing full recovery.
Inducing alternans.#
Using one of the mutant release units, we observe the following dynamics:
By replacing the release unit, we observe the 2:1 alternating behaviour. The release unit in this case releases a different amount of water (its the two holed mutant shown above).
Although it quite easy to carry out these type of experiments, It quickly becomes rather complicated to measure a full set of initial conditions and parameter values. For instance to have a cheap and on the fly way to measure the water flow, or to automate the process of water recycling over many tries becomes pretty complicated if we are running the experiments on a near to zero budget. This makes simulation a great tool to use.
Simulations.#
Interactive simulator: The tunable parameters are H: Red bar position, H2: Blue bar position, Xo: Stimulus angle range (white line), km: Fraction of water released, Is: current strength. The plots displayed are the mass present at the right end of the seesaw \( \Delta m (t) \) as well as the angle temporal evolution \( \Phi(t) \). Once a solution is computed the PLAY button shows the rod oscillation simulation.
Playing around with the integrator, We can always find solutions which resemble the experiments. In the case of the recovery from inactivation, it is notable that the mechanism translates in a straightforward way!I have used this small system to connect introduce concepts in genetics, synthetic biology, differential equations, computational physics and recently in physics informed machine learning to audiences of a very varied nature, and thought It would be nice to have it more or less documented in an organised and interactive way accessible for future use.
This is the first out of two entries in the series, the other one linking chemical reactions and electronic circuits.
References.#
- Are SR Ca content fluctuations or SR refractoriness the key to atrial cardiac alternans?: insights from a human atrial model.- AJP. Heart and Circulatory Physiology.
- Minimal model for calcium alternans due to SR release refractoriness.- Chaos.