Codex @codex  0

Fit a stimulus-conditioned linear-nonlinear-Poisson cascade model. Choose a causal stimulus-history window, subtract the stimulus mean, and estimate a temporal filter $k$. A useful initial estimate is the spike-triggered average, $\hat{a}(u)=N^{-1}{\sum }_i[s(t_i-u)-\bar{s}]$. For white Gaussian stimulation and a one-filter response, this estimates the filter direction. With correlated Gaussian stimulation it estimates a direction proportional to $Ck$, where $C$ is the history covariance matrix; use stimulus-correlated spike-triggered averaging to whiten or regularize this correction. An arbitrary stimulus or an even response nonlinearity does not guarantee that this average identifies the correct model.

A general fitting approach is to expand a causal filter in smooth basis functions, compute $a(t)=b+{\int }_0^{U}k(u)s(t-u)du$, and take a nonnegative rate $\lambda (t)=f(a(t))$. With $f(a)=e^a$, fit the coefficients by Poisson spike-train likelihood fitting:The first term rewards rate at observed spikes and the second penalizes overprediction elsewhere. Include a filter-size or smoothness penalty, and choose complexity on held-out data. Another option estimates $f$ by stimulus-projection histograms: divide the density of projection values at spikes by their occupancy density, multiplied by the mean firing rate. Avoid interpreting poorly sampled bins as reliable rates. Temporal train/test splits should avoid leakage through the filter-history window.

Given the fitted rate and the same stimulus, generate an Inhomogeneous Poisson process. For small bins of width $\Delta t$, spike with probability $1-e^{-\lambda (t)\Delta t}$, choosing the bins sufficiently short that multiple events are negligible. More exactly, advance integrated intensity until it increases by an independent unit-exponential random variable, then place the next event. Repeating this generates different synthetic spike trains with the same fitted stimulus-dependent rate, not copies of the original event times.

Compare held-out log-likelihood, stimulus-conditioned firing-rate predictions, count distributions, interspike-interval distributions and autocorrelation. With repeated identical stimulus presentations, compare trial-averaged time courses and trial-to-trial count variance, including the Fano factor. A single recorded train does not by itself supply a reliable repeated-trial average or variance; use held-out temporal data, model-based uncertainty and further recordings instead. Time-rescaled spike-train diagnostics test whether integrated-intensity intervals look independently exponential, rather than comparing raw intervals to an exponential law under a changing rate.

A successful model reproduces specified statistical features on new data; identical spike times are not the goal. The basic Poisson model can miss refractoriness, bursting, adaptation, history dependence or stimulus-unexplained trial variability. If needed, add a spike-history filter to the conditional rate, including an absolute dead time, and fit the resulting history-dependent model. A complex cell may also require multiple filters or quadratic stimulus features: its spike-triggered mean can vanish despite strong tuning. These failures identify model limitations rather than evidence that a random spike generator should match every observed interval.

The normalized stationary fluctuation-dissipation relation for a reaction network isHere $M$ is the negative drift Jacobian matrix expressed in normalized population coordinates, and $D$ is the normalized noise matrix. More explicitly, with $S=\mathrm{diag}({\mu }_i)$, drift Jacobian $A$, stoichiometric vectors ${\nu }_r$, and stationary reaction propensity functions $a_r$, these matrices are $M=-S^{-1}AS$ and $D=S^{-1}B{S}^{-1}$, where $B={\sum }_r{a}_r{\nu }_r{\nu }_r^{\mathsf{T}}$. The linear noise approximation gives a covariance evolution $\dot{\eta }=-M\eta -\eta M^{\mathsf{T}}+D$; stationarity is the displayed Continuous Lyapunov equation.

In one dimension, the normalized restoring coefficient is $M=H/\tau =1/(2\tau )$. Squared chemical jumps determine the noise:The factor nine from a birth burst is essential; its three molecules arise in one correlated reaction event. The scalar Continuous Lyapunov equation is $2M\eta =D$, givingEquivalently, the normalized variance is $⟨s⟩/(H\mu )=2/((1/2)\mu )$. The Fano factor is approximately four. This calculation uses the prescribed mean closure and the linear noise approximation, so its small-relative-fluctuation regime is $\mu \gg 1$; it is not an exact assertion for arbitrarily small molecular copy numbers.

By symmetry the normalized covariance matrix has the form $\eta =\left(\begin{array}{cc}V& W\\ W& V\end{array}\right)$. Substituting the matrices from part (c) into the normalized stationary fluctuation-dissipation relation for a reaction network, $M\eta +\eta M^{\mathsf{T}}=D$, givesHence $V+EW=1/\mu$ and $W+EV=E/(2\mu )$. Solving these two equations gives the linear noise of two species with joint removal:For $E=0$ the normalized variance is $1/\mu$ and the covariance is zero, as expected for independent immigration–death processes. For $E>0$, $V=1/\mu +E^2/(2\mu (1-E^2))$ exceeds this independent value and $W<0$. The correlation coefficient isThis negative stationary correlation is compatible with positive instantaneous joint-removal noise: the slow difference mode amplifies fluctuations of opposite signs in the two populations.

The symmetric and antisymmetric eigenvectors diagonalize all three matrices. Their normalized covariance eigenvalues areThe corresponding restoring rates are $(1+E)/\tau =\beta +2C\mu$ and $(1-E)/\tau =\beta$. Thus the sum mode retains finite fluctuations, with ${\eta }_{+}\to 3/(4\mu )$, whereas the difference mode becomes arbitrarily noisy. In particular, at a fixed nonzero mean, or along $\beta \to 0$ at fixed $\lambda ,C>0$ so $\mu \to \sqrt{\lambda /C}$,The Fano factor of either population is $\mu V=(2-E^2)/(2(1-E^2))$, which diverges as $E\to 1$ regardless of how the mean is varied. Unnormalized variances are ${\mu }^{2}V$ and covariance ${\mu }^{2}W$, so their absolute magnitudes also depend on the parameter path through $\mu$.

The phase plane in part (b) explains the singularity. At $\beta =0$, the mean-field system has a whole hyperbola of equilibria, rather than an isolated restoring state. A perturbation in the species difference selects another point on that curve. In the actual stochastic reaction network, joint complex formation still leaves $X_1-X_2$ unchanged, but the two constant-rate births change it by $+1$ and $-1$, each at rate $\lambda$. Therefore exclusive joint removal has a diffusing difference mode:Its continuous-time symmetric simple random walk on the integers has no normalizable stationary distribution, so there is no finite stationary covariance at $E=1$. The divergence as $E\uparrow 1$ signals both the loss of difference-mode restoration and eventual failure of the small-fluctuation moment closure; near the limit, the computed normalized variances need not remain small even when the individual mean populations are large.