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 oscillating system 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 below (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 detailed model of atrial cell electrophysiology. (Left): Human atrial cell compartment model schematics, (Right): The RyR2s states Markov model with the four configurations O: open, C: closed, and Inactive states I1 and I2.
The diagram above shows several cell compartments connected through pumps and gates, as well as the connections of several inner and outer domains of the cell (Left panel). And the dynamics responsible for the calcium cycling, where the Ryanodine receptors (RyR2) states can be found in four configurations or states \( (I_1, I_2, O,C) \). The reciprocal of the activation rate \( \tau_r = \frac{1}{k_B} \) corresponds to the recovery from inactivation time.
The translation of the sketched currents 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 is a simplified and elegant version that serves very well to illustrate what we want in this post. The diagram of the currents consideres is shown below:
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 equations used in [2] read as follows:
$$\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 the rhs of (1) is the external stimulus \(I_{CaL}(t)=I_{s}(t)\), Which is provided by the pacemaker cells in the heart or externally in experiments.
A healthy cell will normally release Calcium periodically, driven by the periodic stimulus. If the recovery time of the ion release channels is fast compared to the stimulus period. 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 slow due to a malfunction of the calcium release units, the periodic response alternates between a high/low beat to beat calcium release. 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 as τrec varies. From here we can find the intervals which 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 toys. It is important to take into consideration the system important parameters.

To incorporate the stimulus, we need to add mass to the initial mass attached to the right end of the rod. A way to do this is to use the function:
$$m_2=m_{2o} + \alpha t$$
If the angle \(\Phi \geq \Phi^{\dagger} \). Where \(\Phi^{{\dagger}}\) is a parameter specifying the range of action if the stimulus.
The release 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 remaining in the system after the boundary condition at \( \Phi_2 \) is reached. This parameter is the analogous to the parameter \( k_{im}\). It regulates the recovery time.
Other important parameters are: the threshold mass \(m_1 \) which needs to be surpassed by \(m_2\) to kickstart the motion and the current strength \(\alpha\) specified above.
Experiments.#
The release units can be modified in a modular way, and plugged in to the mechanism, just as we do in some gene regulatory networks. The mutants basically consist of release mechanisms which loss water, faster or slower depending on the release position and number of holes in the recipient.
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 refill 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 dynamics 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.