What should a derivation explain?

A small burst enters a quiet network. With one balance of excitation and inhibition, the burst dies away. With another, activity sustains itself. The goal is to locate the boundary from the network’s equations, before deciding what its avalanche distributions should look like.

The starting ingredients are an equation describing how activity changes and a parameter we can vary. We first find activity levels that remain constant. Then we ask what a small disturbance does near each level. Finally, we locate where that behavior changes and check what the full equation predicts there. The result is a candidate boundary, together with the assumptions under which it was derived.

This tutorial follows that route for continuous-time dynamics. Later, interacting populations introduce oscillations, and approximations introduce a separate question: how accurately do the analyzed equations describe the network? Familiarity with a derivative as a local rate of change is enough to begin. A derivation of a stability boundary is one route into criticality; it is not a universal test for every phase transition.

One activity level, two competing effects

Consider a nonnegative amount of activity. At low activity, recruitment can exceed loss or fall short of it. At higher activity, crowding limits further growth. We call the activity $x$, time $t$, and the adjustable low-activity growth setting $r$. Here these are scaled model quantities, not neuron counts or laboratory seconds.

For a concrete example, take activity 0.1. At growth setting 0.4, recruitment contributes 0.04 while crowding subtracts 0.01: activity initially increases. At growth setting −0.4, both terms reduce activity. The following rule captures those two effects:

We can now draw two connected pictures. The first follows activity as time passes. The second asks, for each possible current activity, how fast does it change? On that second graph, a point above the horizontal zero line means activity increases; a point below it means activity decreases. Its horizontal axis is activity itself.

Start at r = 0.4 and time 0. The gold point is x = 0.1 in the first graph. At the same activity in the second graph, F(0.1, 0.4) = 0.03. That positive height gives the upward slope of the time curve at that instant. Move the time slider to follow this connection. Changing r instead changes the rule governing the entire trajectory.

The trajectories suggest a boundary at r = 0. To derive it rather than read it from a picture, first ask which activity levels can remain constant.

Where can activity remain constant?

The zero-height points on the rate-of-change graph identify activity levels that remain constant. A silent system stays silent under this rule: with no activity, there is nothing to recruit more activity. At r = 0.4, activity 0.4 also stays constant because recruitment and crowding exactly balance. A state that stays constant is a steady state, also called a fixed point. We mark one with a star, $x^*$.

A branch is a family of solutions as the parameter varies. The second branch represents positive activity only when r is positive; a negative algebraic solution lies outside this model’s physical domain. We will call zero the silent steady state, and the positive solution the active steady state. We next examine the active state at r = 0.4 before following what happens as r decreases.

We now need to distinguish a steady state the system returns to from one it moves away from.

Does activity return after a small disturbance?

Finding the critical setting requires distinguishing two behaviors: a small amount of activity fades toward silence, or it grows and sustains itself. The steady-state equation tells us where activity can remain constant. We now need a way to determine whether activity returns after a small displacement. First we will calculate that response at one fixed setting; then we will lower r to find the boundary between returning to silence and sustaining positive activity.

Keep r = 0.4 and consider the active steady state x* = 0.4. Move activity to 0.45. Our rule gives F(0.45, 0.4) = 0.18 − 0.2025 = −0.0225: activity decreases toward 0.4. At 0.35, F(0.35, 0.4) = 0.14 − 0.1225 = +0.0175: activity increases toward 0.4. Because both nearby starting points move back, this steady state is locally attracting. A steady state that amplifies small displacements away from it is unstable.

We want one calculation that captures this local behavior. Measure activity relative to the steady state: the difference is a perturbation, written η (eta), also denoted Δx: the change in activity from x*. At activity 0.45, η = +0.05; at 0.35, η = −0.05.

On the graph of F against x, those two nearby points lie on a downward-sloping curve through (0.4, 0). Its local slope captures the restoring response. We calculate that slope by differentiating F with respect to activity x, keeping r fixed:

To draw the straight line that matches the curve locally, we have a point, (0.4, 0), and a slope, −0.4. This line is the tangent. We know the horizontal displacement η = Δx = x − 0.4 and want to find the corresponding vertical change, Δy. Here the vertical coordinate represents the predicted rate of change of activity. On a straight line, slope is vertical change divided by horizontal change:

Multiply both sides by Δx to solve for Δy. With slope −0.4 and displacement Δx = η = 0.05:

The tangent starts at vertical value zero, so its value after this displacement is 0 + (−0.0200) = −0.0200. This is its predicted rate of change of activity. Calling the tangent’s vertical value L(x), we can write the same construction for any nearby x:

For this displacement η = 0.05, the tangent’s −0.0200 is close to the full equation’s −0.0225. Using this line to approximate F near the steady state is linearization. Move the displacement below to compare both predictions at the same activity.

At r = 0.4, the active steady state is x* = 0.4 and its local slope is −0.4: small displacements return to that positive activity. This establishes that the network can sustain activity at this particular r. Our goal, however, is to derive the critical point: the value of r separating sustained positive activity from eventual decay to silence. Lowering r weakens the growth term rx. Starting from the active state we have just studied, how far can we lower r before sustained positive activity is lost? Locating that boundary will give us the critical point; the slope calculation will help us establish the behavior on either side.

How far can r decrease while activity sustains itself?

The steady-state equation already gave us the positive solution x* = r. Lowering r from 0.4 to 0.2 therefore moves the active steady state from activity 0.4 to 0.2. At r = 0.1, it is 0.1. We name this branch explicitly to distinguish it from silence:

The sustained activity approaches zero as r approaches zero. Below zero, x* = r would be negative, outside our model of nonnegative activity. This suggests r = 0 as the boundary we want to derive: the value of r where the network changes from returning to silence to sustaining positive activity.

What happens to the silent steady state during this change? Substitute zero activity into the original equation, keeping r arbitrary:

Exactly zero activity stays zero for every r. To distinguish the two outcomes, start with the same small positive activity, x = 0.1, and compare different values of r. At r = −0.4 it dies out; at r = 0 it also dies out; at r = 0.4 it approaches the active steady state 0.4. The arrows below show how activity moves along the activity axis.

We need to calculate the r value separating these outcomes. The local slope tells us whether a small positive displacement from zero initially decreases or increases. We call the slope at the silent steady state a(r); r is now the parameter we vary, while the point at which we evaluate the slope stays x = 0.

Differentiate F with respect to x to obtain r − 2x, then substitute x = 0 to evaluate the slope at the silent steady state:

For a small positive x, the tangent predicts a rate of change rx. At r = −0.4 that change is negative: activity fades. At r = +0.4 it is positive: activity grows away from silence. The dividing setting is where this slope passes through zero. Solving a(r) = 0 gives r = 0, the candidate critical setting.

We must also check that the full equation supports this transition. The graphs keep the initial activity at 0.1 while r changes. Compare −0.4, 0, and +0.4: below zero the trajectory approaches silence; above zero it approaches a positive activity. This connects the slope calculation at silence to the behavior of the whole trajectory.

Starting from the same activity x = 0.1, the system settles toward silence when r ≤ 0, and toward a positive active steady state when r > 0. At r = 0, the silent state changes from attracting positive activity to repelling it as r increases, and a positive attracting steady state appears. A parameter change that alters the existence or stability of steady states in this way is called a bifurcation. Here it marks the transition from decay to silence to sustained activity—the critical point we set out to derive.

At what coupling can neural activity sustain itself?

In our first example, r directly controlled growth near silence. For neurons, we want to obtain that growth rate from the processes that activate and deactivate them. Our goal stays the same: find the coupling at which a small amount of activity can sustain itself without continuing external input.

Consider neurons that switch between inactive and active states. Write their active fraction as S. If S = 0.1, ten percent are active and ninety percent are available to activate. Each inactive neuron activates at a rate determined by its input; each active neuron deactivates at rate α. Multiplying each rate by the fraction eligible for that change gives the gain and loss of active fraction.

In this population model, all neurons receive the same input. Existing activity contributes w₀S, where w₀ is the net recurrent coupling. External drive contributes h. The activation-rate function f converts their sum into an activation rate:

As in the earlier example, first establish that the equation can sustain positive activity. Set external drive h = 0 so that recruitment comes entirely from existing activity. For a concrete activation function, call the input s and take f(s) = tanh(s) for positive input, and zero for nonpositive input. This function starts at zero, rises with input, and levels off. Keep the deactivation rate α = 0.1. In these model units, the complete equation is:

Begin at coupling w₀ = 0.2 and initial activity S = 0.1. The time curve rises toward an active steady state near S = 0.498. The adjacent rate curve explains why: below that steady state its height is positive, so activity increases; above it the height is negative, so activity decreases. This is the same restoring behavior we found around x* = 0.4 in the previous model.

Our goal is to locate the critical coupling separating sustained activity from decay to silence. Lower w₀ using the controls below, keeping the initial activity at 0.1. The positive steady state moves toward zero. By coupling 0.1, the initial activity decays toward silence. The graphs suggest a boundary; we will now derive its location.

Previously we varied r; here we vary w₀. Previously x evolved; here S evolves. To locate the boundary analytically, we again calculate the slope of the rate curve at the silent steady state while varying the parameter.

First establish that silence remains a steady state throughout the comparison. With external drive h = 0 and f(0) = 0, substituting S = 0 gives:

Once every neuron is inactive, the zero-drive system cannot restart itself. This makes silence an absorbing state. The boundary is where a small positive activity changes from decreasing to increasing. We therefore need the slope of the full rate curve at S = 0, rather than its value there, which is already zero for every coupling.

To calculate that slope, follow what a small activity does: it supplies recurrent input, that input activates neurons, and some active neurons deactivate. We already know the input is s = w₀S. What activation rate does such a small input produce?

The activation curve below answers that question. Keep w₀ = 0.2 and consider a small active fraction S = 0.05. Its input is s = 0.2 × 0.05 = 0.01, giving f(0.01) = tanh(0.01) ≈ 0.0099997. The straight line through zero with slope 1 predicts 0.01. Near zero input, this tangent closely follows the activation curve.

Call the activation function’s positive-input slope β. In this example β = 1. For a different activation function its value could differ. The tangent’s height is slope × input, so it predicts the activation rate βs:

Now substitute s = w₀S. This turns the tangent prediction βs into βw₀S: activity produces input, and input produces an activation rate.

Multiply this activation rate by the inactive fraction 1 − S to obtain recruitment, then subtract deactivation αS. Expanding makes the term proportional to S visible, which is the tangent contribution we need:

The term proportional to S gives the tangent at silence. Its coefficient is the local slope a, just as r was in the first example. Recruitment contributes βw₀; deactivation subtracts α. We can now solve for the coupling where these balance:

Here w₀c names the candidate critical coupling. With α = 0.1 and β = 1 in compatible model units, it is 0.1. At coupling 0.05, the slope is 0.05 − 0.1 = −0.05, so a small positive activity decreases. At coupling 0.2, it is 0.2 − 0.1 = +0.1, so a small positive activity increases.

The calculated boundary w₀c = 0.1 agrees with the change in the full-equation plots above. Below it, positive activity dies away; above it, the attracting positive steady state supports sustained activity. This is the same sequence as the earlier derivation: establish the active state, follow it toward silence, and locate the boundary using the slope at silence.

We have derived a boundary for one shared activity variable. When excitatory and inhibitory populations have different dynamics, activity in one changes the input to the other. We next need to follow their responses together.

How do we check several activities at once?

We located the neural threshold by asking whether a small activity grows or fades. With interacting populations, that growth can come from a loop: extra activity in the first recruits the second, which recruits the first in return. Can this feedback amplify a disturbance even when each population would restore its own activity in isolation?

To isolate that question, consider two populations near a chosen steady state. Write their departures from that state as η₁ and η₂. Zero displacement means each population is at its steady activity. Give only the first population an initial displacement of 0.02; leave the second at zero.

In this illustrative local model, each population restores its own displacement with rate contribution −2 times that displacement. The other population contributes c times its displacement, where c is the mutual recruitment strength. We will vary c while keeping the restoring rate and initial disturbance fixed. These two positively coupled populations let us examine feedback before adding inhibition.

The arrows below show the loop. At weak coupling, the second population initially rises in response to the first, but both departures eventually fade. At strong coupling, recruitment around the loop wins and both departures grow. The plots show departures from steady activity, so returning to the horizontal zero line means returning to that steady state.

Checking only each population’s own restoring contribution, −2, would give the same answer in both panels. To locate the coupling where the outcome changes, we must include the incoming contribution too. For the first population, add −2η₁ and cη₂; for the second, add cη₁ and −2η₂:

There are four local responses to keep track of. Arrange them with the receiving population’s rate as the row and the changed population as the column:

This matrix is the Jacobian. Each entry is a local slope. For example, increasing η₂ by 0.01 increases the first population’s rate by c × 0.01, giving the top-right entry c. In a full model, we obtain these entries by differentiating each population’s rate with respect to each activity at the chosen steady state.

We now have the local equations, but still need the coupling at which disturbances start growing. In the strong-feedback plot, the two populations increasingly move together. To calculate the growth rate of that coordinated motion, consider equal displacements in both populations. This is one contribution to the original disturbance; the remaining contribution will be accounted for below.

For this contribution, take a displacement of 0.01 in each population. Substitute those two displacements. Because they are equal, each rate becomes the same multiple of 0.01:

At c = 1, both displacements shrink. At c = 3, both grow. The multiplier c − 2 is playing the role of our scalar slope a. The vector (1, 1) specifies the relative displacements—both activities changing together—and the matrix returns a scaled copy of that vector. A vector with this property is an eigenvector; its multiplier is the corresponding eigenvalue.

We must also check disturbances that change the activities differently. For this matrix, the independent direction (1, −1) has eigenvalue −2 − c: one activity rises while the other falls, and that difference fades for c ≥ 0. Every two-entry disturbance can be built from these two directions. Thus both directions fade for 0 ≤ c < 2; the together direction starts growing for c > 2. Solving c − 2 = 0 locates the stability boundary at c = 2.

The feedback plots are therefore separated by c = 2: below it every small disturbance fades; above it the coordinated disturbance grows. This locates a loss of stability of the chosen steady state. The full population equations must determine where the growth eventually leads. So far both populations recruit each other. If the second instead suppresses the first, feedback can make activity swing above and below its steady value. How do we judge whether those swings are returning or growing?

What if activity oscillates while returning?

Suppose a rise in excitation recruits inhibition. Inhibition then pushes excitation down; as excitation falls, inhibition also subsides, allowing excitation to recover. A disturbance can travel around this loop repeatedly. Its sign changes during that motion, so an upward or downward movement by itself no longer tells us whether the system is returning to steady activity.

We need to compare the size of successive swings. The plots below start from the same disturbance as before: η₁ = 0.02 and η₂ = 0. They show the first population’s departure from its steady activity. In the left panel each swing becomes smaller; in the right each becomes larger. Both complete a cycle in the same time.

Separate two visible features: how quickly the swings repeat, and how their size grows or shrinks. The repetition rate stayed fixed between the panels; the change in size switched from decay to growth. To connect that observation to a stability calculation, use an illustrative local feedback rule. The second population suppresses the first with contribution −η₂, while the first recruits the second with contribution +η₁. Each also has its own contribution a times its displacement. We vary a and keep those cross-population influences fixed:

At the initial disturbance (0.02, 0), the second population starts increasing at rate 0.02. Its growing η₂ subtracts from the first population’s rate. Following the two equations produces the plotted swings. Their solution separates the repeating motion from its changing size:

Cosine and sine describe motion around the steady state, completing one turn every 2π time units. The common factor 0.02eᵃᵗ gives its distance from that state. At a = −0.1 this distance shrinks; at a = +0.1 it grows. The zero line in the plots is the same steady state throughout.

What do the eigenvalues report for this matrix? They are λ = a ± i, where i is the imaginary unit, satisfying i² = −1. The real part a gives exactly the growth rate of the distance we just followed. The imaginary parts ±1 describe its rotation. Thus these eigenvalues capture both features visible in the plots.

To locate the stability boundary, vary a while keeping the feedback loop fixed. Negative a gives shrinking oscillations; positive a gives growing oscillations. The dividing value is a = 0. At that value the local linear equations predict a constant distance from the steady state; the full equations determine what the remaining terms do.

For a loop with another oscillation rate, write the eigenvalues as λ = a ± iω. Here ω is the angular frequency: one cycle takes 2π/|ω| time units. The same growth criterion applies. For smooth continuous-time equations, negative real parts for every eigenvalue imply local attraction; a positive real part permits growth. The candidate oscillatory boundary is:

The full equations determine whether growing swings settle into a sustained oscillation. When a pair of complex eigenvalues crosses this boundary under the conditions described below, a nearby family of periodic motions emerges. This is a Hopf bifurcation. It locates an onset of oscillation, whereas the zero-drive scalar example located an onset of sustained positive activity.

Both calculations locate a change in the behavior of the equations we analyze. Before assigning that boundary to a simulated neural network, we need to establish how those equations were obtained.

How does the population equation relate to the neurons?

Our neural derivation began with one active fraction S and a shared input w₀S. A simulation instead keeps track of individual neurons and their events. To connect the two, we must explain why the population’s gain and loss can be represented by that shared-input equation.

If every inactive neuron has the same activation rate, counting inactive neurons gives the activation term (1 − S)f(w₀S + h). When neurons receive different inputs, they can have different activation rates. Replacing those individual inputs by a population average is an additional approximation: the average of their activation rates can differ from the activation rate evaluated at their average input.

A description built from population averages and assumptions about their interactions is called a mean-field description. The assumptions determine which variation between neurons it retains and how it treats fluctuations of the population activity itself. A stability calculation then predicts the boundary of that population description. Comparing with individual-neuron simulations checks how well it predicts the network’s behavior.

The mean-field approximation tutorial develops this reduction further. Here the practical consequence is to carry the population equation’s assumptions into the interpretation of its critical point.

Use the derived boundary as an experimental reference

For the zero-drive neural example, we now have a concrete prediction: with the stated activation and deactivation rules, the population equation changes from decay to silence to sustained activity at w₀ = α/β. We can choose couplings on either side, introduce the same initial activity, and compare the resulting network recordings with that prediction.

A finite stochastic network produces different trajectories across repeated trials. Even when its population equation supports positive steady activity, a finite zero-drive realization can eventually reach extinction. The simulation comparison therefore needs stated population sizes, trial durations, and repeated realizations. These specify how we observe the predicted transition.

Once the reference is established, avalanche analysis and PCA/RG ask what signatures the recordings exhibit around it. For a co-transmission comparison, derive the reference for each architecture, then compare the range of parameters over which the chosen signature appears. This lets us distinguish movement of the theoretical critical point from a broader range of measured signatures around it.

A compact reference for the derivation

For the onset studied here, specify the zero-drive activity equation, establish silence as a steady state, calculate growth near silence, solve for the parameter where it changes sign, and check that the full equation supports sustained activity above that boundary. For interacting populations, the Jacobian and its eigenvalues supply the local growth calculation.

State x; steady state x*

The current activity and an activity level whose net rate of change is zero.

Silent and active steady states

The silent state has zero activity for every r. The active state has activity r and is positive when r > 0.

Control parameter r; dynamics F

The setting varied in the toy equation and its complete net rate of change.

Tangent L(x)

The straight-line approximation to F near the chosen steady state.

Perturbation η; local slope a

A small departure from a steady state and its first-order growth or decay rate.

Jacobian J; eigenvalue λ; eigenvector v

The matrix of local slopes, a direction’s growth rate in continuous time, and the vector specifying that direction.

Coupling c

The adjustable mutual recruitment in the illustrative two-population local equations.

Angular frequency ω; phase θ; amplitude A

The rotation rate, position within a cycle, and size of an oscillation.

Neural fraction S; α, β, w₀, h

Active fraction, deactivation rate, positive-input activation slope, net coupling, and external drive.

Input s; activation-rate function f

The input w₀S + h received by a neuron and the function converting that input into its activation rate.

Conductance g; voltage V; reversal potential E

The factors determining voltage-dependent synaptic current.

Continue with mean-field approximation for the reduction from neurons to populations, or the Widom line for response peaks in driven systems.