The model of (Equation 1) is very similar to a piecewise linear caricature of FitzHugh-Nagumo (FHN) neuron model (Gerstner and Kistler, 2002), also called the McKean's caricature (McKean, 1970; Tonnelier, 2003). Mathematically, the FHN model is given by:

where f(u) is a polynomial of third degree, e.g., f(u) = u − u³/3, and I_ext is the parameter for bifurcation, as opposed to g_s in Equation (1). In the FHN model, one variable (u), possessing cubic nonlinearity, allows regenerative self-excitation via a positive feedback, and the second, a recovery variable (w), possessing linear dynamics, provides a slower negative feedback. It was reasoned in McKean (1970) that the essential features of FHN model are retained in a “caricature” where the cubic non-linearity is replaced by a piecewise linear function f(u). Nullclines of (Equation 2) with a piecewise linear f(u) are shown in Figure 2B in the phase space u × (1 − w). A function f(u) is trivially possible such that it is equal to v_dd − h(i_i, s) in the regions M and I, hence making the u-nullcline similar to the i_i-nullcline in those regions. In the region N, the difference between f(u) and v_dd − h(i_i, s) for any state s does not result in a difference in the direction of system trajectories but only in their velocity, because for small L/C the trajectories are almost parallel to x-axis. Bifurcation in VO₂ neuron is achieved by tuning the load line using a tunable resistance (g_s), or a series transistor (Figure 3A). Figure 3B shows two load line curves corresponding to different gate voltages (v_gs), where one gives rise to spikes while the other results in a resting state.

Moreover, a single dimensional piecewise approximation of the system can be performed using a dimensionality reduction by replacing the movement along the eigenvector parallel to the x-axis with an instantaneous transition from A to A′, or B to B′. This leaves a 1-dimensional subsystem in M and I each along the V-I curves AB′ and BA′. Experiments using VO₂ show that the metallic state conductance g_vm is very high which causes the charging cycle of v_o to be almost instantaneous (Figure 4) and resembles a spike of a biological neuron. As such, the spiking statistics can be studied by modeling just the discharge cycle of v_o. The inductance being negligible can be effectively removed and only the capacitance is needed for modeling the 1D subsystem of insulating state (Figure 6A) making v_i = v_dd − v_o.

Analytical expressions for the FPT of OU process (with μ = 0) for a constant boundary were derived using the Laplace transform method in Ricciardi and Sato (1988). Reproducing some of its results, let the first passage time for the system (Equation 3), with μ = 0, which starts at x(0) = x₀ and hits a boundary S, be denoted by the random variable t_f(S, x₀), and its mth moment by τ_m(S, x₀). Also, let tf˜(S,x0) be the FPT for another OU process with μ = 0, θ = 1, and σ = 2, and τm~(S,x0) be its mth moment. Then time and space scaling for the OU process imply that

where α=2θσ2. The first two moments for the base case OU process τ1~ and τ2~ are given by

where ϕ_k(z) can be written as an infinite sum

with ρ(n, k) being a function of the digamma function (Ricciardi and Sato, 1988).

We extend this framework for calculating the FPT statistics with a fluctuating boundary S as follows. Let the IMT threshold be represented by the random variable v_h. For the VO₂ based IMT neuron, the 1D subsystem in the insulating phase can be converted in the form of Equation(3) with μ = 0 by translating the origin to the fixed point. If this transformation is T then x = Tv_i = T(v_dd − v_o), S = Tv_h, and x_o = Tv_l. The start and end points are B′ and A, respectively in Figure 2. v_h is assumed constant during a spike, and across spikes the distribution of v_h is v_h ~ D, where D is either Gaussian, or subGaussian. For subGaussian distributions we use the Exponential Power family EP[κ], κ being the shape factor. Let the interspike interval of IMT neuron be denoted by the marginal random variable timt(D,vl). Then t_imt is related to t_f in Equation (4), given common parameters θ and σ, as follows:

The moments of t_imt can be calculated as:

where α=2θσ2. If D is Gaussian or EP[κ] distribution and αT is an affine transformation, then αTv_h also has a Gaussian or EP[κ] distribution.

First moment of t_imt is calculated using Equations (5) and (7) as

The expansion for ϕ_k(z) in Equation(6) can be used to calculate 𝔼_vh[ϕ_k(αTv_h)] using the moments of αTv_h as follows

Figure 6B shows firing rate (1/𝔼[timt(D,vl)]) as a function of v_gs for various σ_t values and for three distributions of threshold fluctuations. The calculations approximate the experimental observations well for all three v_h distributions, the closest being EP[3] with σ_t = 4.

For higher moments, higher order terms are encountered. For example, in case of the second moment, using Equations(5) and (7), we obtain

with a higher order term ϕ1(αTvh)2. In the case of the third moment we obtain ϕ₁(αTv_h)ϕ₂(αTv_h). As each ϕ_k term is an infinite sum, we construct a cauchy product expansion for the higher order term using the infinite sum expansions of the constituent ϕ_ks and then distribute the expectation over addition. For example, if the ϕ_k expansions of ϕ₁(z) and ϕ₂(z) are (∑a_i) and (∑b_i), respectively, then the cauchy product expansion of ϕ₁(z)ϕ₂(z) can be calculated as ∑c_i, where c_i is a function of a_1…i and b_1…i, and the expectation 𝔼[ϕ₁(z)ϕ₂(z)] = ∑𝔼[c_i]. Since c_i is a polynomial in z, 𝔼[c_i] can be calculated using the moments of z.

If μ_imt and σ_imt are the mean and standard deviation of interspike intervals t_imt, the coefficient of variation (σ_imt/μ_imt) varies with the relative proportion of the thermal and the threshold induced noise. Figure 7 shows σ_imt/μ_imt (calculated using parameters matched with our VO₂ experiments) plotted against σ_t for various kinds of v_h distributions fitted to experimental observations. σ_imt/μ_imt as observed in our VO₂ experiments is about an order of magnitude more than what would be calculated with only thermal noise using such a neuron, and hence, threshold noise contributes significant stochasticity to the spiking behavior. As the IMT neuron is setup such that the stable point is close to the IMT transition point (Figure 3B), low σ_t results in high and diverging σ_imt/μ_imt for any distribution of threshold noise, and σ_imt/μ_imt reduces with increasing σ_t for the range shown. For a Normally distributed v_h the variance diverges for σ_t ≲ 8, but for Exponential Power (EP) distributions with lighter tails, the variance converges for smaller values of σ_t. Statistical measurements on experimental data, as indicated in Figure 7, provide measures of σ_imt/μ_imt (dotted line) and σ_t (shaded region). We note that EP distributions provide a better approximation of the stochastic nature of experimentally demonstrated VO₂ neurons as the range of σ_t is estimated to be <5.