Griffiths phase and long-range correlations in a biologically motivated visual cortex model

Activity in the brain propagates as waves of firing neurons, namely avalanches. These waves’ size and duration distributions have been experimentally shown to display a stable power-law profile, long-range correlations and 1/f b power spectrum in vivo and in vitro. We study an avalanching biologically motivated model of mammals visual cortex and find an extended critical-like region – a Griffiths phase – characterized by divergent susceptibility and zero order parameter. This phase lies close to the expected experimental value of the excitatory postsynaptic potential in the cortex suggesting that critical be-havior may be found in the visual system. Avalanches are not perfectly power-law distributed, but it is possible to collapse the distributions and define a cutoff avalanche size that diverges as the network size is increased inside the critical region. The avalanches present long-range correlations and 1/f b power spectrum, matching experiments. The phase transition is analytically determined by a mean-field approximation.

. Elements of the V1 model. A: A realization of the L = 5 network. Circles are photoreceptors (Input), neurons (internal layers) and axonal compartments (Output). Layers are linearized. B: Architecture of the network. C: Spatial organization of the network of N = 4L 2 neurons. The columnar structure is highlighted in red. There is a column of size N c = 4l 2 = 196 neurons centered on each neuron of the network. D: Compartmental scheme of neurons. The probability, P(a k ), of choosing a presynaptic axonal compartment, a k , of any neuron is exponential such that most of the synapses start from the end of the axon (left). The probability, P(d m ), of choosing a dendritic postsynaptic compartment, d m , is Gaussian with mean 50 and standard deviation 10, so that most of the synapses lay in the middle of the dendrite (right).
P G (x, y; x j , y j , σ d ), centered in the (x j , y j ) coordinates of the presynaptic neuron with standard deviation σ d = 3 in each direction. The connections from Input layer to LGN are created using the same procedure.
A columnar structure emerges in the network as a consequence of picking neighbors in the limited region of l 2 neurons. V1 displays such structure as seen in experiments. 4 Therefore, each column has approximately N c = 4l 2 = 196 neurons, see From \ To LGN VI IVCβ II/III Output  Input  1  ----LGN  --500  --VI  --1100  350  -IVCβ  -600  -700  -II/III  ----100  Table S1. The quantity of attempted synapses per presynaptic element from a specific layer (rows) to another layer (columns). Presynaptic element may be a photoreceptor (from Input layer) or a neuron (from internal layers). There are 100 photoreceptors in the Input layer for each neuron of the LGN. There are 100 axon compartments of each neuron in the II/III layer connecting to V2. These values are predefined in the beginning of the simulation, but are not exactly achieve in practice because of free boundary conditions. To LGN VI IVCβ II/III Output  Input  98  ----LGN  --488  --VI  --1073  293  -IVCβ  -585  -683  -II/III  ----1  Table S2. The quantity of synapses per postsynaptic element that a specific layer (column) recieves from another layer (row). Postsynaptic elements are all neurons. There are 100 axonal compartments in the output layer for each neuron and each of them forms a synapse with V2. The averages were performed when assembling a L = 99 network and are similar to the values from Table S1. exponential probability P E (k) = (10/4) exp (10k/4), of the presynaptic cell, j, and one dendritic compartment, d The probability of choosing a dendritic compartment is taken as Gaussian because consecutive events of connecting to any two compartments should be uncorrelated and most of the synapses should arrive around the center of the dendrites. 5 The probability of choosing an axonal compartment is exponential because axons should transport signals for as long as possible, thus most of the outward synapses must come from the end of the axon. [6][7][8][9] The neurons of this model could be mapped into simple celular automata, like the one studied by Kinouchi & Copelli. 10 One should then adjust the probability of exciting a neighbor (in the referred model) as proportional to the fraction of the signal that would reach the soma in our model (derived in Section 3 of this Supplementary Material). The spike would have to be delayed in order to mimic the average amount of compartments in the neuron and the amount of refractory states would have to be adjusted accordingly in order to avoid self-sustained activity in the interlayer loops. Then, it is just a matter of building the structure of the network according to the rules presented in this section.

Macroscopic measurements
The activity A(t) of the network at time t is just the sum of all neuronal firings at each time step, where δ i, j is the Kronecker delta. The temporal profile of A(t) is presented in Figure S2. The density of activated neurons is simply where T is the processing time of the network. Notice that ρ is already normalized as each neuron can only fire once. The avalanche size, s(n), is the sum of network activity, where t n is such that A(t n ) = 0 ∀ n. Figure S2 (bottom) shows a detail of A(t) for E = 1.2 mV. In this figure, we can see that A(t) is zero for some time steps in between peaks. Such intervals are used to separate consecutive avalanches (i.e. consecutive peaks) in our data. The processing time is T = max n {t n } ; the largest avalanche is M = max n {s(n)} (the average is taken across several trials for each E) and the ratio m = M/N is the fractional largest avalanche. The avalanches autocorrelation function, C(t ′ ), and power spectrum, S( f ), are defined as: The characteristic time, τ, of the autocorrelation function is defined by the exponential cutoff fit of The detrended fluctuation analysis (DFA) is also a measure of long-range correlations for time series. [11][12][13] We apply this technique to the activity time series, A(t), in order to obtain its fluctuations. First we obtain the average where T is the processing time. Then, we perform a cumulative sum over the detrended time series: for t = 1, · · · , T . We split the total time interval T in many boxes, each of size ∆t. We then calculate the local trend ofĀ(t) inside each window ∆t simply by fitting a linear function toĀ(t) inside ∆t. Let y ∆t (t) be the y-coordinate of this fit, we may calculate the fluctuations by subtractingĀ(t) from each local trend inside each ∆t window:

Average signal in soma
Consider the structure described in Section 1 of this Supplementary Material for the dendrites of the neurons and the dynamics described in the main text. If a single synapse arrives at the compartment d n (0) of some neuron at instant t = 0, then after k − n time steps the signal in dendrite d k (t = k − n) is d k = λ k−n E with average d k = λ k−n E (we omit time for simplicity). The synapses arrive in any dendritic compartment, d n , with Gaussian distribution P G (n; 50, 10) with mean d (c) = 50 and standard deviation σ d = 10. Then, the average λ k−n is given by: where the integral is continuously evaluated within infinity limits because more than 99% of the Gaussian distribution lies within the bounds of the dendrite. The integrand in Eq. 10 may be manipulated in order to still get a Gaussian function multiplied by a constant: whereP G (n) is a normalized Gaussian with displaced mean. Notice that r k (λ ) ≡ λ k−n = d k /E is the fraction of the signal that reaches compartment k. If the synapse (or external signal) arrived just at compartment d (c) (i.e. σ d = 0), then r k (λ ) = λ k−d (c) would be a simple exponential decay with characteristic time −1/ ln(λ ), as experimentally expected. 5 In experiments the arriving signal compartment should be farther than 100 units from the soma. In our model the signal that reaches the soma is the signal from k = 100. The average signal in soma is therefore d 100 = r 100 E.

Maximum Likelihood Estimation of Power Laws
We performed a Maximum Likelihood test in our avalanche distributions power-law fit. The columnar structure of the model is responsible for the presence of a bump in the shape of the obtained avalanche distributions. To take this into account we used a censored Maximum Likelihood estimator which uses only the data that exhibits a power law in order to calculate the best fit exponent,α. The formula follows the description for discrete power laws made by: 14 whereα is the expected value of the exponent of the best fit power law, s ≥ s min ≥ 1 and the Hurwitz zeta function is defined by where

5/9
We applied this test to the complementary cumulative distributions for some E values and L = 99. For E = 1.10 mV, E = 1.19 mV (critical point) and E = 13.0 mV we obtained best fit slopes with valuesα = 1.46,α = 1.34 andα = 1.50, respectively. For E = 1.88 mV we divided the distribution in two parts in order to avoid the bump around avalanches of size N c ≈ 200. For the first part (with s min = 1, s max = 190) we obtained a slope with valueα 1 = 1.69 and for the second one (considering s min = 600, s max = 1500) we obtained the valueα 2 = 1.57. All these values are fairly close to the values obtained by the avalanche cumulative distribution collapse and F (s) = c 1 + c 2 s −α+1 fit presented in the main text.

Finite-Size Scaling of the fractional largest avalanche
The fractional largest avalanche (or cluster), m, is sometimes regarded as an order parameter to study the Griffiths phase. [15][16][17][18] In Figure S3 we show m as function of E for increasing L and the FSS of m for several E.
Inside the critical region we expect that m ∼ L −D ′ as it is a possible order parameter. Indeed, our fit for 1.11 ≤ E ≤ 1.19 mV (inside the Griffiths phase) gives D ′ = 1.0(1). This value allows us to calculate the scaling of the largest avalanche (or cluster), M, because m ≡ M/N and N = 4L 2 . Therefore,

Finite-Size Scaling of autocorrelation characteristic time
A widely used criterion to determine whether the system presents a Griffiths phase is the FSS of the autocorrelation characteristic time, τ, of the avalanche (or cluster) time series. 15-17 A system with a single critical point and no Griffiths phase is expected to have τ scaling as where L is system size and z is known as the dynamical exponent and is generally non-universal. 19 However, inside a Griffiths phase, τ behaves as a stretched exponential: 15-17 where c is just a constant and D is expected to be the exponent of the largest avalanche defining a typical dimensionality of the system. Figure S4 shows the same simulation data for τ(L) fitted by either the stretched exponential of Equation (19) (panel A) or the power law of Equation (18) (panel B). We ignored τ(L = 20) for fitting the stretched exponential. Notice that both panels present plausible fits, although the stretched exponential fit yields D = 1.1(1) which matches our expectations from other independent calculations of D (from the avalanches cutoff fit, Z, presented in the main text and from the largest fractional avalanche fit, m, presented in the previous section).

Processing time
The network processing time has been discussed in the main text and is again presented in Figure S5A for completeness. The FSS of T ∼ L µ is fitted in Figure S5B and gives µ = 0.80 (1). The variance of T is defined as: Inside the Griffiths region we verified that Var(T ) behaves as: We present Var(T ) as function of E in Figure S5C and its FSS fit, Figure S5D, yields µ ′ = 3.5 (2). There is a striking similarity between Var(T ) and χ ρ ≡ NVar(ρ)/ ρ , specially concerning the exponent values (µ ′ and γ/ν ⊥ ).
For E > E c = 1.19 mV and every L the processing time asymptotic behavior may be fitted by a crossover of two PL (see Figure S5A): where ω = 3.1 (2). The term E −1/ω dominates for E ≫ E c . This equation is empirically obtained and has no underlying principle. However, it was useful to estimate the simulation time. The understanding of the functional form in equation (22) will be considered in future works. Notice that the dependence on L is explicitly given in the Results section of the main text.

Small size effects on the calculations
Ignoring small sized networks simulation data points for fitting is generally acceptable due to finite-size effects (i.e. boundary effects) that are enhanced in small systems. In our case, each column in a network with L = 20 includes a fraction of more than 12% of all neurons in the layers. This is significantly different from fractions of other columns in networks with different L: 3% for L = 40 and less than 1% for L = 80, 99. Avalanches thus are able to grow larger and faster for networks with L = 20 than expected for other L with fixed E.

7/9
Also, initial conditions are a flash of a fixed size square of size 30 × 30 in the Input layer of the model, corresponding to an average activation of 3 × 3 = 9 neurons in the LGN layer. Thus, an average fraction of 9/400 ≈ 2 × 10 −2 of the LGN neurons are initially activated for L = 20. For larger L we have average fractions of 10 −3 or less neurons of the LGN which are initially activated.
The case L = 20 would have the largest avalanche and the activated density slightly larger than expected. On the other hand, the variance of such order parameters has slightly lower values. L = 20 would also present processing time slightly lower than expected because the columns of such network are larger fraction of its layers. If this fact is taken into account and we ignore L = 20 in the FSS analysis of ρ, χ ρ and Var(T ), the exponents β /ν ⊥ , γ/ν ⊥ and µ ′ significantly change from the values presented in the main text and in the previous section to β /ν ⊥ = 0.40(2), γ/ν ⊥ = 2.6(2) and µ ′ = 2.9 (2). The scaling exponent of T , remains practically unchanged µ = 0.81(1) with or without L = 20. However, only the universality class of the model is changed. The phase transition is still a continuous phase transition, such that there is still a Griffiths phase in the very same region of EPSP, 1.11 ≤ E ≤ 1.19 mV.
The best approach to tackle this problem would be to estimate the deviations caused by the small-size effects on L = 20 network and correct its corresponding data points, instead of deliberately ignoring them for fitting. Nevertheless, we do not have a systematic way of accomplishing this, so we are not yet able to definitively determine the model's universality class.  (21), yielding µ ′ = 3.5 (2). Notice the similarity between the FSS of χ ρ (which is the variance of ρ) and Var(T ). (▽) inactive phase, ( ) critical phase, ( ) and (△) active phases. Vertical bars are standard deviation, light gray is the Griffiths phase, dark gray is the weakly active phase and dotted lines are only guides to the eyes.