Bayesian inference for SDE models: a case study for an excitable stochastic-dynamical model

Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer

Abstract: The FitzHugh-Nagumo (FHN) model is a versatile yet simple model for the evolution of excitable or oscillatory physical processes. It is a well-known two-dimensional simplification of the Hodgkin-Huxley model of spike generation in squid giant axons used in neuroscience. The dynamics of the process are governed by a two-dimensional hypo-elliptic Stochastic Differential Equation (SDE) of which only the first component is assumed to be observed discretely in time. Inferring parameters appearing in the SDE is a nontrivial problem. We illustrate how statistical inference can be done within the Bayesian setup using the Backward Filtering Forward Guiding algorithm as detailed in [1]. The complex interaction between the estimated parameters and the dynamics of the model is explored through an interactive interface is phase-space . Due to the generality of the statistical procedure, Julia-implementation and graphical interface, the program can be extended or adapted to other interesting models and applications.


Recurring natural events such as West Antarctic ice shelf collapses, the firing of neurons or the rapid population expansions of defoliating insects typically occur neither in a completely periodic nor completely irregular fashion. The cause is the interplay between the noise coming from outside factors interacting with the system in a random or non-predictable way and the rhythm created by the delayed interaction between excitatory and inhibitory components of the system. Stochastic models are mathematical objects which couple rhythm and noise and can express a wide variety of possible recurrence patterns. Statistical inference for these models cannot be done using standard regression techniques, as these neglect serial correlation and nonlinearities induced by the model. Here, we demonstrate how statistical inference can be performed for a well known model with excitatory and inhibitory components: the FitzHugh-Nagumo (FHN) model ([2], [3]). The procedure takes as input discrete data observed over time and returns samples from the posterior distribution of the parameters that govern the dynamics of the model. We provide an interactive visualisation to assess the effect of each of the parameters on the paths of the diffusion.

The FitzHugh-Nagumo model

The FitzHugh-Nagumo (FHN) model is a versatile yet simple model for the evolution of excitable or oscillatory physical processes. It is a well-known two-dimensional simplification of the Hodgkin-Huxley model of spike generation in squid giant axons used in neuroscience. Mathematically, the model is expressed through a 2-dimensional stochastic system of the form

formula not implemented

where inline_formula not implemented is a scalar Wiener process and inline_formula not implemented are the parameters of the model. Hence inline_formula not implemented is a diffusion process where the first and second components model the membrane potential and a recovery variable respectively. We will assume only the membrane potential inline_formula not implemented is observed at times inline_formula not implemented . Besides intrinsic noise in the equation for the recovery variables, we additionally assume extrinsic Gaussian noise on the discrete-time observations. Hence, we assume observations inline_formula not implemented where

formula not implemented

The SDE for inline_formula not implemented together with the observation scheme for inline_formula not implemented defines a continuous-discrete state-space model. We will consider two related problems:

  • parameter estimation;

  • reconstruction of the latent paths of the diffusion (the smoothing problem).

For the latter problem, the highly nonlinear nature of the drift-term in the FHN-model prevents applications of the Rauch-Tung-Striebel smoother. Instead, we will illustrate how recent research on conditional simulation of diffusions (as described in [1]) provides an elegant solution to the problem. The corresponding framework is generally applicable to SDE-models, but here we solely focus on the FHN-model here for illustration purposes. The statistical method has been implemented in Julia and packaged in Bridge and BridgeSDEInference.

We first ensure all required packages for the analysis ahead are called.

]up; add StaticArrays Distributions Bridge BridgeSDEInference#nextjournal
using Bridge
using StaticArrays
using BridgeSDEInference
using BridgeSDEInference: EulerMaruyamaBounded
using Random, LinearAlgebra, Distributions
const State = SArray{Tuple{2},T,1,2} where {T};

The qualitative behaviour of the solution to the FHN-SDE can be rather different, depending on the chosen combination of parameter values. Without the Wiener noise, sets of parameter values yielding excitatory, oscillatory or "Van der Pol" behaviour are detailed for example in [2] (pages 399-400). Below we illustrate excitatory and oscillatory behaviour for the FHN-SDE considered here.

using Plots
param = :regular # regular parametrization of the trajectories (no transformation of the coordinate processes)
x0 = State(-0.5, -0.6)
# time grid 
dt = 1/1000
T = 20.0
tt = 0.0:dt:T
# Exitatory behaviour (stable fixed point)
ε = 0.1 ; s =0.5 ; γ =1.5 ; β = 1.4 ; σ =0.3 ;
P = FitzhughDiffusion(param, ε, s, γ, β, σ);
X1, _ = simulate_segment(0.0, x0, P, tt);
skip = 10
p1_1 = plot(first.(X1.yy)[1:skip:end], last.(X1.yy)[1:skip:end], leg = false)
x1 = -1.2:0.01:1.2
x2 = [-1.2, 0.3]
x3 = [-0.9, 0.6]
plot!(x1, [x1 .- x1.^3 .+ s])
plot!(x2, γ.*x2 .+ β)
p1_2 = Plots.plot([1:skip:end],  first.(X1.yy)[1:skip:end], leg = false)
Plots.plot!(p1_2,[1:skip:end],  last.(X1.yy)[1:skip:end])
p1 = plot(p1_1, p1_2, layout = (1,2))
# Oscillatory behaviour (unstable fixed point)
β = 0.6
P = FitzhughDiffusion(param, ε, s, γ, β, σ);
X2, _ = simulate_segment(0.0, x0, P, tt); 
p2_1 = plot(first.(X2.yy)[1:skip:end], last.(X2.yy)[1:skip:end], leg = false)
plot!(x1, [x1 .- x1.^3 .+ s])
plot!(x3, γ.*x3 .+ β)
p2_2 = Plots.plot([1:skip:end],  first.(X2.yy)[1:skip:end], leg = false)
Plots.plot!(p2_2,[1:skip:end],  last.(X2.yy)[1:skip:end])
p2 = plot(p2_1, p2_2, layout = (1,2))
p3 = plot(p1, p2, layout = (2,1))
png(p3, "/results/out.png")

Figure: Phase planes (left) and trajectories (right) of two simulations up to time inline_formula not implemented. In both simulations inline_formula not implemented are the same for both the simulations while inline_formula not implemented (top panels) and inline_formula not implemented (bottom panels). The blue and red curves in the right-hand figures correspond to the inline_formula not implemented and inline_formula not implemented components respectively. The green and red curves in the left-hand figures are nullclines.

Data generation

Here we simulate data from which we will extract discrete time observations. To this end, we first simulate a path of the FHN-diffusion on a fine time grid. The parameters are chosen such that the model is sensitive to changes in the parameters (near bifurcation):

# starting point under :regular parametrisation
x0 = State(-0.5, -0.6)
ε = 0.1 ; s =-0.8 ; γ =1.5 ; β = 0.0 ; σ =0.3 ; 
P = FitzhughDiffusion(param, ε, s, γ, β, σ);
# time grid 
dt = 1/50000
T = 20.0
tt = 0.0:dt:T
X, _ = simulate_segment(0.0, x0, P, tt); 

From the simulated path, we retain only a subset of it, reflecting the realistic assumption that the process is partially observed. In particular, as specified above, we assume to observe only the first coordinate inline_formula not implemented at discrete time points with small extrinsic noise:

# subsampling
num_obs = 100
skip = div(length(tt), num_obs)
Σ = [10^(-4)]
L = [1.0 0.0]
obs = (time =[1:skip:end],
  values = [Float64(rand(MvNormal(L*x, Σ))[1]) for x in X.yy[1:skip:end]]);

Whereas we use simulated data here to verify performance of the proposed methods, in any real world application the data are given. These should then replace obs by importing the data as NamedTuple{(:time, :values)}. The observations and the (latent) simulated paths can be visualised as follows (in practise, only the marked observations are given and the paths are unobserved):

skip = 50
Plots.plot([1:skip:end],  first.(X.yy)[1:skip:end], label = "X")
Plots.plot!([1:skip:end],  last.(X.yy)[1:skip:end], label = "Y")
p = Plots.scatter!(obs.time, obs.values, markersize=1.5, label = "observations")
png(p, "/results/out.png")

Figure: Simulated trajectory of the FitzHugh-Nagumo model and observed points (dotted in green). A small perturbation of inline_formula not implemented can result in a large, non-linear excursion of inline_formula not implemented, a so called spike.

Statistical inference

To infer the parameters and reconstruct (smooth) the latent trajectories of the diffusion, we consider a a Bayesian approach using the Backward Filtering Forward Guiding algorithm as detailed in [1]. This is a Markov Chain Monte Carlo algorithm to solve the smoothing problem, i.e. to reconstruct the latent trajectories of the diffusion. This is combined with a data-augmentation approach where the algorithm iteratively updates the latent diffusion path and parameters, conditional on the observations.

In the FHN-model, simultaneous estimation of all the parameters of the model is not possible due to identifiability issues. Therefore in the following we fix inline_formula not implemented to the value used during data generation.

Backward Filtering Forward Guiding

We now sketch the main idea behind this algorithm, referring the reader to [1] for full details. The aim of this algorithm is to do smoothing, i.e. reconstruct the continuous-time path based on the incomplete discrete-time observations. It turns out that the diffusion, conditioned on all information provided by the data, satisfies a stochastic differential equation of the form

formula not implemented


formula not implemented

Hence, the SDE for the conditioned process is obtained from the SDE of the unconditioned process by superimposing a guiding term which depends on inline_formula not implemented. Whereas the precise form of inline_formula not implemented is somewhat daunting, what's most important here is that it is determined by the observations and the intractable transition densities inline_formula not implemented of the diffusion inline_formula not implemented. As inline_formula not implemented is intractable, we replace it with a proxy which is tractable. As in [1], we replace inline_formula not implemented by the transition densities inline_formula not implemented of an auxiliary process inline_formula not implemented, the choice of which we detail below. Let inline_formula not implemented be defined in terms of inline_formula not implemented , just as inline_formula not implemented is defined in terms of inline_formula not implemented. Then, instead of forward simulating inline_formula not implemented, we forward simulate the process inline_formula not implemented defined by

formula not implemented

The discrepancy between inline_formula not implementedand inline_formula not implemented can be adjusted for by the likelihood ratio between their induced laws. So all that needs to be provided is the dynamics of the auxiliary process inline_formula not implemented. As in [1], we take it to be a linear (Gauss-Markov) diffusion. That is

formula not implemented

To reflect the dynamics in the FHN-model in between any two observations inline_formula not implemented, we take the linear diffusion corresponding to the linearisation of the original model at the point inline_formula not implemented. Therefore we set

formula not implemented

Configuration of the algorithm

In the configuration of the algorithm all parameters are initialised. Here we start with a random perturbation of the parameter values used for generating the data:

# Take the real β, as it is fixed. 
θ_init = [ε, s, γ, β, σ].*(1 .+ (2*rand(5) .- 1).*[0.2, 0.2, 0.2, 0.0, 0.2]);

We define the target process FitzhughDiffusion and on each segment between spanning consecutive observations times an auxiliary process FitzhughDiffusionAux (inline_formula not implemented ):

P_trgt = FitzhughDiffusion(param, θ_init...)
P_aux = [FitzhughDiffusionAux(param, θ_init..., t₀, u, T, v) for (t₀,T,u,v)
        in zip(obs.time[1:end-1], obs.time[2:end], obs.values[1:end-1], obs.values[2:end])]
# Container
model_setup = DiffusionSetup(P_trgt, P_aux, PartObs());

We set the observation scheme and the imputation grid dt of the simulated latent path and specify a Gaussian prior on the starting point x0 with mean zero and covariance matrix equal to inline_formula not implemented.

# Observation scheme
L = @SMatrix [1. 0.]     
Σ = @SMatrix [10^(-4)]      
set_observations!(model_setup, [L for _ in P_aux], [Σ for _ in P_aux],
                  obs.values, obs.time)
# Imputation grid
dt = 1/200
set_imputation_grid!(model_setup, dt)
# Prio distribution on (X_0, Y_0) 
set_x0_prior!(model_setup, GsnStartingPt(x0, @SMatrix [.1 0; 0 .1]), x0);

To compute the term inline_formula not implemented, certain systems of ordinary differential equations (ODEs) need to be solved numerically. For that purpose we specify the ODE-solver.

initialise!(eltype(x0), model_setup, Vern7(), false, NoChangePt(100))
# Further setting 
set_auxiliary!(model_setup; skip_for_save=2, adaptive_prop=NoAdaptation());

The final step is to set prior distributions on the parameters and specify the transition kernels for the Metropolis-Hastings update steps for each parameter in the model. Here, we simply take improper (uniform) priors and random walk updates on the parameters with the inline_formula not implemented-distribution. In case the parameter is known to be strictly positive, the update is performed on the log of the parameter. The latter is specified by setting the second argument in UniformRandomWalk to true.

mcmc_setup = MCMCSetup(
      Imputation(NoBlocking(), 0.975, Vern7()),
      ParamUpdate(MetropolisHastingsUpdt(), 1, θ_init,
                  UniformRandomWalk(0.5, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 1))
      ParamUpdate(MetropolisHastingsUpdt(), 2, θ_init,
                  UniformRandomWalk(0.5, false), ImproperPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 2))
      ParamUpdate(MetropolisHastingsUpdt(), 3, θ_init,
                  UniformRandomWalk(0.5, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 3))
      ParamUpdate(MetropolisHastingsUpdt(), 5, θ_init,
                  UniformRandomWalk(0.5, true), ImproperPosPrior(),
                  UpdtAuxiliary(Vern7(), check_if_recompute_ODEs(P_aux, 5))
schedule = MCMCSchedule(10^4, [[1, 2, 3, 4, 5]],
                        (save=10^3, verbose=10^3, warm_up=100,
                         readjust=(x->x%100==0), fuse=(x->false)));

Running the algorithm

The central function for the inferential procedure is mcmc. The function takes as input the configuration we set up above and, while running, prints some useful information regarding the MCMC algorithm and the (default) adaptive scheme which tunes the step sizes of the Markov chains for each parameter.

out = mcmc(mcmc_setup, schedule, model_setup);

Evaluation of the posterior

The output out contains iterates of the posterior distribution of parameter and paths. First we inspect the parameter iterates in the Markov chain:

θ1_names = ["eps" "s" "gamma" "sigma"]
chain2 = out[2].θ_chain[1:10:end]
p = plot(hcat( [i[1] for i in chain2 ], 
  		[i[2] for i in chain2], [i[3] for i in chain2], 
  		[i[5] for i in chain2]), label = θ1_names,layout = (4,1))
png(p, "/results/out.png")

Figure: Trace-plots of MCMC chain for parameters inline_formula not implemented.

From these trace-plots it appears that after approximately inline_formula not implemented iterations the chains have reached their stationary region. We also visualise the smoothed latent paths

ss = 5
p1 = scatter(obs.time, obs.values, color = :green, markersize=1.5)
Plots.plot!(p1,[1:ss:end],  first.(X.yy)[1:ss:end], alpha=0.5, lw = 1.8, color=:green, label = "Xtrue")
Plots.plot!(p1,[1:ss:end],  last.(X.yy)[1:ss:end], alpha=0.5, lw = 1.3, color=:green, label = "Ytrue")
for j in 1:length(out[1].paths)
		Plots.plot!(p1, out[1].time[1:ss:end],  [x[1] for x in out[1].paths[j]][1:ss:end], 
    	color = :blue, alpha = 0.1, leg = false)
  	Plots.plot!(p1, out[1].time[1:ss:end],  [x[2] for x in out[1].paths[j]][1:ss:end], 
    	color = :red, alpha = 0.1, leg = false )
png(p1, "/results/out.png")

Figure: Samples of the inferred latent path, displayed together with the true path in green (i.e. the path from which the discrete time data were extracted).

From that we can extract the 90% credible sets and the median value of the marginal posterior densities of the parameters and visualise the marginal densities of the parameters through boxplots. For that we consider the initial inline_formula not implemented iterations as burnin samples.

# Save median and 90% credible sets of the empirical posterior distribution of the parameters as a dataframe
using CSV
using DataFrames
burnin = 10^4
θ_names = ("eps","s","gamma","beta","sigma")
θ = [ε, s, γ, β, σ]
df = DataFrame(par = String[], true_val = Float64[], median = Float64[], Q1 = Float64[], Q99 = Float64[])
for j in 1:5
  push!(df, [θ_names[j], θ[j], median([i[j] for i in out[2].θ_chain[burnin:end]]),
      quantile([i[j] for i in out[2].θ_chain[burnin:end]], 0.01),
      quantile([i[j] for i in out[2].θ_chain[burnin:end]], 0.99)])
CSV.write("results/posterior.csv", df);
df[!,2:5] = round.(df[!,2:5]; digits = 3)
5 items
using StatsPlots
p = []
for j in [1, 2,3,5]
  push!(p, StatsPlots.histogram([i[j] for i in out[2].θ_chain[burnin:end]], label = θ_names[j]))
png(plot(p[1], p[2], p[3], p[4], layout = (2,2)), "/results/out.png")

Figure: Boxplots of posterior samples of the parameters inline_formula not implemented.

Clearly, there is considerable more uncertainty on inline_formula not implemented and inline_formula not implemented compared to inline_formula not implemented and inline_formula not implemented.

Interface for interpreting the fitted model

(Run from here to get the animation without redoing the statistical experiment)

Here, we are mostly interested in visualising and exploring the interaction between the parameters and the dynamical patterns created by the model and finally reconnect the inferential results with the natural phenomenon. This can be achieved by a cointegration of the Julia package Makie.jl and JavaScript which enable an interactive interface on a web based platform. We fix first some values which determine the initial position of the sliders and their limits:

# Point estimates
using CSV
thetas = round.(
)[!, :median]; digits = 2)
# Needed for sliders as well
thetalims = [(0.0,1.0), (-5.0,5.0), (0.0,10.0), (0.0,10.0), (0.0,1.0)]
truethetas = [0.1, -0.8, 1.5, 0.0, 0.3];
#thetas = truethetas # should true parameters or point estimates marked in the app?

Next, we draw the histograms relative to the posterior distribution obtained from the inferential procedure explained above. The images will lay on the slider and facilitate exploring the behaviour of the model as a function of the parameters:

# Histogram plots to be used as slider background
using Plots
using Plots.PlotMeasures
pl = Vector{Any}(undef, 5)
for j in [1, 2, 3, 4, 5] 
  	pl[j] = Plots.scatter([thetas[j]], [2.5],
        showaxis=false, grid=false, legend=false, border=false, size=(230,100),
        markersize=2.5, color=:black, xlim=thetalims[j], ylim=(0,100), left_margin=-9mm, margin=-5mm, background=RGB(0.95,0.95,0.95), background_color_outside=:white)
  	if @isdefined out
    qs = quantile([i[j] for i in out[2].θ_chain[burnin:end]], 0.0:0.1:1.0)
        Plots.plot!(pl[j], Bridge.piecewise(SamplePath(qs, 0.1*(thetalims[j][2]-thetalims[j][1]) ./diff(qs)))...,
            linewidth=0, color=nothing, fillrange = 0, fillalpha = 0.5,
      	fillcolor = :red)
		savefig(pl[j], "results/slider$j.png");

Finally we are ready to construct the interactive interface which consist of showing the effect on the evolution of the particles when changing the parameters of the model. Each particle represents a possible value of the quantity of interest and the latent variable as point in the 2-dimensional phase plane inline_formula not implemented. The particles move according to the dynamics of the model through this plane. The values of the parameters can be changed through the sliders. Changing the parameters affects the dynamics of the particles, making spikes or excursions more/less likely.

] add Hyperscript Markdown JSServe WGLMakie#master AbstractPlotting Colors;
using JSServe, WGLMakie, AbstractPlotting
using JSServe: JSServe.DOM, @js_str, onjs
using Colors
using Random
using WGLMakie: scatter, scatter!
# fallback values if statistical analysis is not run
if !@isdefined thetas
   thetas = [0.1, -0.8, 1.5, 0.0, 0.3]
if !@isdefined thetalims
   thetalims = [(0.0,1.0), (-5.0,5.0), (0.0,10.0), (0.0,10.0), (0.0,1.0)]
using Hyperscript, Markdown
using JSServe, Observables
using JSServe: Application, Session, evaljs, linkjs, div, Asset
using JSServe: @js_str, onjs, Button, TextField, Slider, JSString, Dependency, with_session
using JSServe.DOM
# load histogram images to use as slider background 
sliderbg = [
styles = map(1:5) do i
  return css("#slider$(i)",
    paddingTop = "10px",
    height = "50px",
    backgroundSize = "115px 50px",
    backgroundRepeat = "no-repeat",
    backgroundPosition = "center center",
    backgroundImage = sliderbg[i])
function dom_handler(session, request)
  	# fetch initial parameter (initial slider settings)
    eps = thetas[1];
    s = thetas[2];
    gamma = thetas[3];
    beta = thetas[4];
    si = thetas[5];
    rebirth = 0.001; # how often a particle is "reborn" at random position
  	sl = 101 # slider sub-divisions
    # slider and field for sigma
    slider5 = JSServe.Slider(range(thetalims[5]..., length=sl), si)
    nrs5 = JSServe.NumberInput(si)
    linkjs(session, slider5.value, nrs5.value)
    # slider and field for beta
    slider4 = JSServe.Slider(range(thetalims[4]..., length=sl), beta)
    nrs4 = JSServe.NumberInput(beta)
    linkjs(session, slider4.value, nrs4.value)
    # slider and field for gamma
    slider3 = JSServe.Slider(range(thetalims[3]..., length=sl), gamma)
    nrs3 = JSServe.NumberInput(gamma)
    linkjs(session, slider3.value, nrs3.value)
    # slider and field for s
    slider2 = JSServe.Slider(range(thetalims[2]..., length=sl), s)
    nrs2 = JSServe.NumberInput(s)
    linkjs(session, slider2.value, nrs2.value)
    # slider and field for eps
    slider1 = JSServe.Slider(range(thetalims[1]..., length=sl), eps)
    nrs1 = JSServe.NumberInput(eps)
    linkjs(session, slider1.value, nrs1.value)
    # slider and field for rebirth
    slider6 = JSServe.Slider(0.0:0.0001:0.005, rebirth)
    nrs6 = JSServe.NumberInput(rebirth)
    linkjs(session, slider6.value, nrs6.value)
    # init
    R = (1.5, 3.0) # plot area
    R1, R2 = R
    limits = FRect(-R[1], -R[2], 2R[1], 2R[2])
    n = 400 # no of particles
    K = 150 # display K past positions of particle fading out
    dt = 0.0005 # time step
    sqrtdt = sqrt(dt)
		particlecss =
    	css("input[type='number']", width = "135px"), 
 			css("input[type='range']", width="135px")
    ms1 = 0.02 # markersize particles
    ms2 = 0.02 # markersize isokline
  	# plot particles, initially at random positions
    scene = WGLMakie.scatter(repeat(R1*(2rand(n) .- 1), outer=K), repeat(R2*(2rand(n) .- 1),outer=K), color = fill((:white,0f0), n*K),
        backgroundcolor = RGB{Float32}(0.04, 0.11, 0.22), markersize = ms1,
        glowwidth = 0.005, glowcolor = :white,
        resolution=(500,500), limits = limits,
  	# style plot
    axis = scene[Axis]
    axis[:grid, :linewidth] =  (0.3, 0.3)
    axis[:grid, :linecolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.3),RGBA{Float32}(0.5, 0.7, 1.0, 0.3))
    axis[:names][:textsize] = (0.0,0.0)
    axis[:ticks, :textcolor] = (RGBA{Float32}(0.5, 0.7, 1.0, 0.5),RGBA{Float32}(0.5, 0.7, 1.0, 0.5))
    splot = scene[end]
  	# plot isoklines
    WGLMakie.scatter!(scene, -R1:0.01:R1, (-R1:0.01:R1) .- (-R1:0.01:R1).^3 .+ s, color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2)
    kplot1 = scene[end]
    WGLMakie.scatter!(scene, -R1:0.01:R1, gamma*(-R1:0.01:R1) .+ beta , color = RGBA{Float32}(0.5, 0.7, 1.0, 0.8), markersize=ms2)
    kplot2 = scene[end]
  	# set up threejs scene
    three, canvas = WGLMakie.three_display(session, scene)
    js_scene = WGLMakie.to_jsscene(three, scene)
    mesh = js_scene.getObjectByName(string(objectid(splot)))
    mesh1 = js_scene.getObjectByName(string(objectid(kplot1)))
    mesh2 = js_scene.getObjectByName(string(objectid(kplot2)))
    # init javascript
    evaljs(session,  js"""
        iter = 1; // iteration number
    		// fetch parameters
        eps = $(eps);
        s = $(s);
        gamma = $(gamma);
        beta = $(beta);
        si = $(si);
        R1 = $(R1);
        R2 = $(R2);
    	  rebirth = $(rebirth);
    		// update functions for isoklines
        updateklinebeta = function (value){
            beta = value;
            var mesh = $(mesh2);
            var positions = mesh.geometry.attributes.offset.array;
            for ( var i = 0, l = positions.length; i < l; i += 2 ) {
                        positions[i+1] = beta + positions[i]*gamma;
            mesh.geometry.attributes.offset.needsUpdate = true;
            //mesh.geometry.attributes.color.needsUpdate = true;
        updateklinegamma = function (value){
            gamma = value;
            var mesh = $(mesh2);
            var positions = mesh.geometry.attributes.offset.array;
            for ( var i = 0, l = positions.length; i < l; i += 2 ) {
                        positions[i+1] = beta + positions[i]*gamma;
            mesh.geometry.attributes.offset.needsUpdate = true;
            //mesh.geometry.attributes.color.needsUpdate = true;
        updateklines = function (value){
            s = value;
            var mesh = $(mesh1);
            var positions = mesh.geometry.attributes.offset.array;
            for ( var i = 0, l = positions.length; i < l; i += 2 ) {
                        positions[i+1] = positions[i] - positions[i]*positions[i]*positions[i] + s;
            mesh.geometry.attributes.offset.needsUpdate = true;
            //mesh.geometry.attributes.color.needsUpdate = true;
    		// move particles every x milliseconds
            function (){
                function randn_bm() {
                    var u = 0, v = 0;
                    while(u === 0) u = Math.random(); //Converting [0,1) to (0,1)
                    while(v === 0) v = Math.random();
                    return Math.sqrt( -2.0 * Math.log( u ) ) * Math.cos( 2.0 * Math.PI * v );
                var mu = 0.2;
                var mesh = $(mesh);
                var K = $(K);
                var n = $(n);
                var dt = $(dt);
                var sqrtdt = $(sqrtdt);
                k = iter%K;
                var positions = mesh.geometry.attributes.offset.array;
                var color = mesh.geometry.attributes.color.array;
                for ( var i = 0; i < n; i++ ) {
                    inew = k*2*n + 2*i;
                    iold = ((K + k - 1)%K)*2*n + 2*i;
                    positions[inew] = positions[iold] + dt/eps*((1 - positions[iold]*positions[iold])*positions[iold] - positions[iold+1] + s); // x
                    positions[inew+1] = positions[iold+1] + dt*(-positions[iold+1] + gamma*positions[iold] + beta) + si*sqrtdt*randn_bm();
                    color[k*4*n + 4*i] = 1.0;
                    color[k*4*n + 4*i + 1] = 1.0;
                    color[k*4*n + 4*i + 2] = 1.0;
                    color[k*4*n + 4*i + 3] = 1.0;
                    if (Math.random() < rebirth)
                        positions[inew] = (2*Math.random()-1)*R1;
                        positions[inew+1] = (2*Math.random()-1)*R2;
                for ( var k = 0; k < K; k++ ) {
                    for ( var i = 0; i < n; i++ ) {
                        color[k*4*n + 4*i + 3] = 0.98*color[k*4*n + 4*i + 3];
                mesh.geometry.attributes.color.needsUpdate = true;
                mesh.geometry.attributes.offset.needsUpdate = true;
        , 15); // milliseconds
  	# react on slider movements
    onjs(session, slider1.value, js"""function (value){
        eps = value;
    onjs(session, slider2.value, js"""function (value){
      onjs(session, slider3.value, js"""function (value){
    onjs(session, slider4.value, js"""function (value){
    onjs(session, slider5.value, js"""function (value){
        si = value;
  	onjs(session, slider6.value, js"""function (value){
        rebirth = value;
  	# arrange canvas and sliders as html elements
    dom = DOM.div(, particlecss, DOM.p(canvas), DOM.p("Parameters"), DOM.table("ε ($eps)"),"s ($s)"),"γ ($gamma)"),"β ($beta)"),"σ ($si)")),, id="slider1"), DOM.div(nrs1)),, id="slider2"), DOM.div(nrs2)),, id="slider3"), DOM.div(nrs3)),, id="slider4"), DOM.div(nrs4)),, id="slider5"), DOM.div(nrs5))
       ),"rebirth", DOM.div(slider6, id="slider6"), DOM.div(nrs6)),
# attach handler to current session
JSServe.with_session() do session, request
    dom_handler(session, request)
Shift+Enter to run

Animation: Particles moving according to the inferred model on the phase plane inline_formula not implemented . The sliders are showing the marginal posterior and are initially set to the median value of the inferred parameters (black dot).


In this article we have considered the FHN-model and shown how parameters and latent paths can be reconstructed from partial observations. Nextjournal offers a flexible framework for reproducible research and sophisticated visualisation. Due to the generality of the statistical procedure, we believe extensions to other interesting models and applications are feasible. This project is embedded in a larger research project on statistical methods for stochastic dynamical systems that team members have been working on over the past 5 to 10 years.


  • [1] M. Mider, M. Schauer and F.H. van der Meulen (2020): Continuous-discrete smoothing of diffusions,

  • [2] P. Sanz-Leon, S.A. Knock, A. Spiegler and V.K. Jirsa (2015): Mathematical framework for large-scale brain network modeling in The Virtual Brain, NeuroImage 111, 385-430.

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


This work has been supported by NextJournal with the grant "Scholarship for Explorable research". We are grateful to NextJournal for their continuous support, for the great computational resources made available to us and for the warm hospitality offered during our visit to NextJournal's headquarters in Berlin.

  • Project Exploring and Statistically Learning an Excitable Stochastic-Dynamical Model. Scholarship for Explorable Research,, 2019-2020.

Runtimes (1)