Coupled FitzHigh - Nagumo System

During the 1950s Hodgkin and Huxley came up with a set of differential equations which mimic the physics of neural activity - based on experiments with the giant squid axon.

In essence, a neuron doesn’t fire until a threshold input stimulus is reached and, after firing, neurons exhibit a relaxation time, which prevents them from firing again for a short time.

A simplified version of the (four dimensional) Hodgkin-Huxley model is the (two dimensional) Fitzhugh -Nagumo [1] oscillator:

In the above, x is the neuron membrane voltage, a, b, c and k are constants and y is the 'recovery variable' for the neuron membrane potential. The membrane voltage (x) has a cubic non-linearity.

Fitzhugh-Nagumo oscillators can be coupled together - Sprott considers the following [2] coupling of two Fitzhugh-Nagumo oscillators - the stimulus for the first is some constant (k) times the membrane voltage of the second (x₂); the stimulus for the second is some constant (k) times the membrane voltage of the first (x₁):

Note: there appears to be a typographical error in Sprott's equation (6.10) - the third equation should (I think) be as I've written above - i.e. y₂ rather than y₁ appearing on its right hand side.

Here's my patching for the coupled equations:

Richard FitzHugh with analog computer, ca. 1960,
from [3].

Analog computer, ca. 2015. Patched for
coupled FitzHugh - Nagumo system.

The following plots were obtained using the parameter values / initial conditions given by Sprott [2] (i.e. a = 0.4, b = 0.5, c = 0.01 and k = -1). The plots show the first ten minutes of the trajectories - using the 0.2 µF integration capacitors throughout. The red plot shows y₂ vs. x₂; the black plot shows y₁ vs. x₁.

10 minutes' worth of trajectory. Green dots mark start / initial condition
point in each case (red pen took time to get going!).

With the computer in FAST mode I captured the following (y₁ vs. x₁)...just before my Tektronix 465 blew up(!).

Attractor for the coupled FitzHugh - Nagumo system,
as per [2], figure 6.7. (1 V/div both axes.)

It's interesting to watch the trajectory unfold with the computer running in slow time (0.2 µF integration capacitors = 0.2 second time constant), on the plotter. The pen whizzes around the perimeter from the initial condition position, settling down at (almost) a point (towards bottom left - resting state - bright splodge on oscilloscope display); the trajectory then slowly spirals outwards (initially very slowly) from this resting point, eventually gaining speed (and eventually reaching the firing threshold / stimulus) and going around the perimeter again (quickly), and so on.

Sometimes the fast part of the trajectory follows the outer path / orbit, followed immediately by a second slightly inner fast path / orbit, before returning to the 'resting spiral' (equilibrium resting potential). The fast moving outer parts of the trajectory corresponds to the neuron firing; the spiral part corresponds to the equilibrium / resting part. I guess the second fast orbit which is, on occasion, observed corresponds to a relative refractory (see below) part of the trajectory.

Reference [3] gives a good depiction of the trajectory and its elements, reproduced below. The neuron fires (active part of trajectory); immediately followed by a part in which it cannot be stimulated, no matter how great a stimulus is applied (absolute refractory), followed by a 'relative refractory' part (in which a second firing is possible - but only given enough stimulus) - followed by the resting phase.

Phase portrait and physiological state diagram of FitzHugh-Nagumo model, taken from [3].

References

[1] FitzHugh R. (1961) Impulses and physiological states in theoretical models of nerve membrane. Biophysical J. 1:445-466.

[2] Sprott J.C. Elegant Chaos: Algebraically Simple Chaotic Flows, pp. 123 - 125.

[3] http://www.scholarpedia.org/article/FitzHugh-Nagumo_model

100 weeks

Almost exactly 100 weeks since I started, I've finished. Never again!

Baltic birch plywood side panels. Recessed
handles by Mentor.

Rear of the computer. Two 120 mm fans at top.

Funny - doesn't look like a USB connector
- because it isn't - this allows a second
computer (or, more realistically, part thereof)
to be connected...

Hour counter, IEC mains connector and fuse.

Conjecture: A tilted circular billiard table is chaotic

Here we have a point-like ball moving inside a (soft, elastic) circular boundary i.e. a circular billiard table. Interestingly, Charles L. Dodgson, (Lewis Carroll of Alice's Adventures in Wonderland fame), once published a set of rules for a two-player game of circular billiards [1]. However, this is predated somewhat by Alhazen / Ptolemy [2].

This post follows the same approach as the previous post, but with the V-wedge replaced with a full circle of radius R:

q is the distance from a point (x,y) to the circular surface, measured along the normal to the surface.

Hence if a point-like ball is at a point (x,y) outside the circular perimeter, there is an elastic restoring force (acting towards the centre of the circle) given by kq, where k is the spring constant of the wall of the circle. It is easy enough to work out q, and thus the accelerations along x and y:

where m is the mass of the ball.

Acceleration due to gravity is directed purely along the negative y-axis. We assume that the radius of the circle is small compared to the radius of the earth (i.e. g is constant).

The computer's squarers yield the square of the input voltage divided by ten - both quantities Q and q end up divided by the square root of ten; this comes out in the wash of the divider, however. The only point to remember is that the voltage corresponding the the radius R must also be divided by √10 also (before it gets subtracted from Q/√10 to yield q/√10). Here's the patching:

No gravity case

We simply set g = 0.

Now the ball doesn't move (there is no initial force). However if we set the initial velocity to something non-zero the point-like ball moves across to the circular wall, and rebounds. This is most interesting perhaps if we start the ball from someplace other than the origin - otherwise the ball simply oscillates back and forth through the origin. The following plots show the trajectory for the ball starting at (1,1) with an initial velocity of 3 m/s along the positive x-axis. -R√10 is set to -0.632 V, corresponding to a billard table radius R of 2 (metres).

In all cases, the trajectory was plotted for eight minutes (integrators have 1 second time constants). Different coloured plots are for different wall spring constants: (a) k = 0.2 N/m (black line), (b) k = 1 N/m (red line), (c) k = 10 N/m (these spring constants assume a 1 kg mass).

The results indicate that is indeed not possible for the trajectory of a moving ball within a circular billiard table to cover the whole table completely. Unless of course...

Gravity Case

...we tilt the table. In other words we have a non-zero gravitational force (here along the negative y-axis).

We start with the ball at (1,1) as before, but with zero initial velocity. Acceleration due to gravity (g)  is 10 m/s². The following result is for k = 10 N/m (eight minutes of trajectory).

A lesser tilt also gives extensive coverage of the (lower half) of the table...here g = 2 m/s² (about an 11° slope). Everything else as before:

An even lesser tilt (g = 0.5 m/s² - about a 2.8° slope on planet Earth) still yields good coverage - albeit at a slower pace:

Finally, g = 0.1 m/s² (about a 0.57° slope on planet Earth - and surely hardly noticeable in a darkened pool hall)  yields the result below. This result is suspect to some degree however, since the ball passes higher than its starting height - suggests some issue with the computation - e.g. an offset null somewhere (e.g. squarer outputs), in effect inputting energy into the system:

References

[1] Dodgson, Charles Lutwidge ('Lewis Carroll'). Circular Billiards for two players. 1890.

[2] Alhazen's problem corresponds to finding the point on the edge of a circular billiard table at which a cue ball at a given point must be aimed in order to bounce off the edge of the table and strike another ball at a second given point, named after the Arab scholar Abū ʿAlī al-Ḥasan ibn al-Ḥasan ibn al-Haytham (c. 965 AD – c. 1040 AD), although it was first formulated by Ptolemy even earlier, in 150 AD.

Gravitational Wedge Billiard

Here we have a V-shaped barrier in which a point mass (ball) bounces about, under the influence of a constant vertically downward force (e.g. gravity). There seemingly being nothing new in chaos, this is (perhaps well) known as the gravitational wedge billiard [1].

We consider a right angle wedge as shown below:

As in the previous post, I've assumed that if the ball is below slope there is an outward (towards positive y-axis) restoring force along the normal to the slope. As before, this force is given by Hooke's law. In the previous post a value for q (the normal distance from ball to slope) was worked out, applicable to the right hand side (i.e. x > 0). Extending the thing to both sides of the wedge gives the result:

and, as before m is the mass of the ball and k the spring constant. Hence the above gives the motion of the ball as it hits, descends into, and eventually bounces off the slope. The bounce is elastic. The bigger the value of k, the less spongy the slope is. Here's the patching:

Each of the following results shows the first four minutes' worth of the trajectory (1 second time constant on the integrators), with the ball starting at (0,6). Acceleration due to gravity (acting vertically downwards) is 10 m/s². The quoted spring constants correspond to a 1 kg mass. Note that, as spring constant k is reduced, the ball sinks further below the x axis, at least for the first few bounces. Although, in theory, one would expect the ball to simply bounce up and down along the y-axis forever, in practice it deviates to one side or the other after two or three bounces, enough to flick it on to one of the sides of the wedge, high up enough to lead to further deviations and ensuing chaos.

k = 200 N/m

k = 20 N/m

k = 10 N/m

An interesting observation: the 'effective' wedge angle reduces as the spring constant reduces: the more spongy the slope, the smaller the wedge angle. What's also interesting / unexpected (perhaps) is that the sides of the triangular wedge trajectory space remains sharply defined.

Finally, if, rather than dropping the ball straight down the wedge's vertical axis, we launch it from part way up the left hand side, towards the right hand side, we achieve a more stable trajectory / orbit. For example, with the initial velocities along x and y as +8 m/s and +6 m/s respectively, and starting at (-4,4) we get:

This looks something like a precessing(?) period two orbit.

References

[1] Classical and quantum chaos of the wedge billiard. I. Classical mechanics, T. Szeredi and D. A. Goodings, Physical review. E, November 1993; 48(5).

Bouncing Ball II

Here, the ball hits a (45 degree) slope.

If the (point-like) ball is at some position (x,y) below the slope, as shown below, then I am assuming an elastic restoring force acting along the normal to the slope, and that this force is proportional to the distance q, according to Hooke's law:

Hence, if we are below the slope (i.e. x > y), we have (where m is mass of ball, k the spring constant of the slope and g the acceleration due to gravity):

And if the ball is above the slope, then we just have gravity acting (along the negative y-axis) on the ball. This is easy enough to patch up:

The right hand side summing amplifier (10x gain) has both a negative (i.e. inverting) and positive (i.e. non-inverting output, both used here. The switch closes every time the ball goes below the surface - thus applying the restoring acceleration. The left integrators also have gains of 10x, and hence k/(2m) = 100 and if we assume a mass of 1 kg, we thus have a spring constant of 200 N/m.

2 V/cm both axes - ball starts at (4,6)

2 V/cm both axes - ball starts at (5,6). Second bounce visible.

The upper trace shows the ball dropping vertically until it hits the slope and then bouncing of the slope. The lower trace shows what happens when the ball is dropped closer to the slope (i.e. further to the right) - the ball has less kinetic energy when it hits the slope the first time and thus makes a lesser bounce. The ball comes down to hit the slope a second time (at about x = -3). In all cases the ball sinks into the slope somewhat - due to the finite and slightly squidgy (200 N/m) spring constant which was assumed.

Bouncing Ball I

I like this sort of high school physics. Made me think what happens when the ball hits the surface (at y = 0). I've included a spring component to give an elastic bounce. When the ball is above the y = 0 plane, the only force acting is gravity and so, in this region,  the acceleration y'' equals g, where g is acceleration due to gravity (e.g. -10 m/s² , minus since gravity is acting downwards). When the ball is below the plane, there is an additional elastic restoring force (acting upwards), equal to -ky (remembering y is negative below the plane!), so that y'' equals g - ky/m in this region.  

A damping force is readily included (representing a simple mechanical viscous damper (dashpot)), equal to -dy', where d is the damping coefficient (newton-seconds per meter), and y' the velocity. The initial height is specified by setting the initial condition of the second integrator (e.g. 5 volts for 5 m).

Initial height = 5 m; k = 100 N/m; d = 0.5 Ns/m.
Vertical axis 2 V per division, horizontal axis 0.2 seconds per division.