Activity propagation with Hebbian learning

TL;DR

This paper explores how Hebbian learning influences infection and activity spreading models, revealing complex behaviors like phase transitions, immune regions, and Griffiths effects driven by local reinforcement mechanisms.

Contribution

It introduces a novel extension of the contact process incorporating Hebbian learning, demonstrating diverse emergent phenomena and global effects from local incentives.

Findings

Negative reinforcement promotes activity spreading and immune regions.

Positive reinforcement can cause Griffiths effects and non-universal scaling.

Emergent behaviors include phase transitions and complex activity patterns.

Abstract

We investigate the impact of Hebbian learning on the contact process, a paradigmatic model for infection spreading, which has been also proposed as a simple model to capture the dynamics of inter-regional brain activity propagation as well as population spreading. Each of these contexts calls for an extension of the contact process with local learning. We introduce Hebbian learning as a positive or negative reinforcement of the activation rate between a pair of sites after each successful activation event. Learning can happen either in both directions motivated by social distancing (mutual learning model), or in only one of the directions motivated by brain and population dynamics (source or target learning models). Hebbian learning leads to a rich class of emergent behavior, where local incentives can lead to the opposite global effects. In general, positive reinforcement (increasing activation rates) leads to a loss of the active phase, while negative reinforcement (reducing activation rates) can turn the inactive phase into a globally active phase. In two dimensions and above, the effect of negative reinforcement is twofold: it promotes the spreading of activity, but at the same time gives rise to the appearance of effectively immune regions, entailing the emergence of two distinct critical points. Positive reinforcement can lead to Griffiths effects with non-universal power-law scaling, through the formation of random loops of activity, a manifestation of the ``ant mill" phenomenon.

Tables

Table 1. Table 1 : Summary of the studied learning models. The precise value of the reinforcement parameter r r is often irrelevant, apart from whether it corresponds to negative ( r < 1 r<1 ) or positive ( r > 1 r>1 ) reinforcement. Gray background color indicates the lack of a phase transition. For r < 1 r<1 , the behavior can depend on the initial condition being a fully active lattice (full) or a single active site (spreading). We use the following abbreviations: DyP: dynamical percolation, GP: Griffiths phase. LDP stands for “learning” DP, where certain aspects are governed by the DP universality class, see text. (*) stands for emergent quenched directional configuration(s) in the r → 0 r\to 0 limit, see text.

	1D		2D
		$r>1$	$r\to 0$	$r\to \infty$
source			full: LDP	GP
	LDP	directional lock	spreading: DyP	(exp. localization)
mutual	full: LDP		full: LDP
	spreading: ballistic	localization	spreading: DyP	localization
target	ballistic front	XXXXX localization	full: LDP spreading:	localization localization
	(*)		DyP (*)

Table 2. Table 2 : Critical scaling in one dimension. Critical activation rate and critical exponents at various negative reinforcement rates ( r < 1 r<1 ). Results on standard directed percolation (DP) are shown up to four digits after the decimal.

	$r$	${\lambda }_c$	$\alpha$	$\delta$	$\eta$	$z$
DP [35, 36]	1.0	3.2979	0.1595	0.1595	0.3137	1.5808
source	0.9	3.512(1)	0.15(1)	0.05(1)	0.45(1)	1.54(3)
	0.8	3.721(1)	0.14(1)	0.04(1)	0.46(2)	1.54(2)
	0.5	4.425(2)	0.15(1)	0.03(2)	0.48(2)	1.53(4)
	0	8.13(1)	0.16(1)	0.04(2)	0.45(2)	1.55(3)
mutual	0.9	3.269(1)	0.29(7)	0.00(1)	0.73(4)	1.005(3)
	0.8	3.248(1)	0.30(8)	0.00(1)	0.75(3)	1.002(4)
	0.5	3.203(1)	0.25(7)	0.00(1)	0.79(4)	0.998(7)
	0	3.193(1)	0.17(2)	0.00(1)	0.83(3)	1.00(2)

Table 3. Table 3 : Critical scaling in 2D. Critical activation rate and critical exponents at negative reinforcement rate r → 0 r\to 0 , in comparison to the directed percolation (DP) and dynamical percolation (DyP) universality classes.

	${\lambda }_1$ (spread)	${\lambda }_2$ (full)	$\alpha$	$\delta$	$\eta$	$z$
DP [37]	1.64874(2)	1.64874(2)	0.4526(7)	0.4526(7)	0.233(6)	1.757(8)
DyP [38, 39]	4.66571(3)	4.66571(3)	0.092	0.092	0.586	1.1293
source	1.6768(1)	1.695(1)	0.44(1)	0.07(7)	0.65(3)	1.12(2)
mutual	1.3990(2)	1.4675(1)	0.44(4)	0.07(4)	0.71(7)	1.08(3)
target	1.5970(5)	1.8360(5)	0.55(4)	0.099(4)	0.57(4)	1.13(2)

Equations

$$N$$

$$P$$

$$R$$

$$\rho \sim t^{-\alpha }$$

$$\kappa =\frac{1-\delta }{2}=0.4203$$

$${\dot{p}}_n=-\mu p_n-\lambda p_{n}q+r_{n-1}{p}_{n-1}q+l_{n+1}{p}_{n+1}q+\lambda p{q}_n.$$

$$\dot{p}=-\mu p+\lambda p(1-p),$$

$${\dot{q}}_n=\mu p_n-\lambda p{q}_n.$$

$$q_n=\frac{\mu }{\lambda }\frac{p_n}{p}$$

$$0=-\lambda p_n+r_{n-1}{p}_{n-1}+l_{n+1}{p}_{n+1}.$$

$$p_n\sim e^{-U_n}\simeq e^{-n^2\frac{\ln (1/r)}{2}}.$$

$$\tau \sim (n^{\ast }{)}^2\sim \frac{1}{r-1}$$

$${\dot{p}}_n=-\mu p_n-\lambda p_{n}q+r_{n-1}{p}_{n-1}q+l_{n+1}{p}_{n+1}q+$$

$$+\overline{rp}p{q}_{n+1}+\overline{lp}p{q}_{n-1},$$

$$0=-2\lambda p_n+\left(r_{n-1}+\frac{\lambda }{2}\right)p_{n-1}+\left(l_{n+1}+\frac{\lambda }{2}\right)p_{n+1}.$$

$$p_n\sim e^{-U_n}\simeq e^{-n\ln 3}.$$

$${\dot{p}}_n=-\mu p_n+\overline{rp}p{q}_{n+1}+\overline{lp}p{q}_{n-1},$$

Full text

Introduction

Biological and social systems show learning and adaptation across a broad range of scales [1, 2, 3]. Yet, understanding the impact of local learning mechanisms remains mostly unexplored even in simple models. Here, we focus on a paradigmatic model of activity (or infection) spreading, the contact process (CP) [4] and explore the impact of local learning on its dynamics, motivated by Hebbian learning [5]. As introduced by T. E. Harris, the CP is a variant of the susceptible-infected-susceptible (SIS) model that is equivalent to the standard SIS model on regular lattices [4]. The CP is a continuous-time Markov process that acts on a network of sites with each site either in an “active” or “inactive” state. Other common terminology for these states are infected/healthy and occupied/unoccupied. Each site can be active multiple times, like how a person can get sick many times in the case of the common cold or COVID, or how a neuron or brain region can be active on multiple occasions.

In the case of infection spreading, it is natural to explore what happens if individuals learn from the dynamics and lower the rate at which they interact with those individuals that infected them, a form of negative reinforcement, or anti-Hebbian learning [6]. For example, a reinforcement effect similar to social distancing [7] would result in a mutual learning mechanism, where the two sites involved would reduce the time spent together.

In the context of brain dynamics, each site can be viewed as a small region of the brain (on the order of a mm3) that can be either active or inactive and each edge can be viewed as a connection between the brain regions [8]. As chemical synapses in the brain are directed, this motivates implementing a directed learning mechanism, either impacting the connection starting at the source or at the target. Positive reinforcement between brain regions aligns with the observation that brain regions that fire together will wire together more strongly [9, 10, 11].

As a further motivation, the CP has been also applied to the spreading of populations [12]. For example, parasites spreading on a background of hosts can induce local adaptations through a triggered immune response and resistance resulting from the host-parasite “arms race” [13, 14, 15, 16], leading to a form of negative reinforcement. Positive reinforcement could also happen along ecological corridors that are strips of favorable conditions (e.g. land and water) connecting fragmented habitats [17]. Once established, these corridors can be preferably used later by the population.

Both negative and positive reinforcement leads to heterogeneities in the system if started from a homogeneous network, a form of kinetic disorder. Here, we address the question of how local learning changes the large-scale emergent behavior of the system, compared to the standard CP. The standard CP has two phases, an active phase and an inactive phase, governed by the activation rate $\lambda$. In the active phase, the infection survives in an infinite system, whereas, in the inactive phase, the infection dies out and never returns—the system reaches an absorbing state. This is known as an absorbing phase transition [18, 19]. The transition in the standard CP serves as a representative of the robust directed percolation universality class (DP) [20, 21], which is believed to contain any absorbing phase transition with short-range interactions, a scalar order parameter, a single absorbing state, and no extra conservation laws or symmetries [22, 23]. Outside of the CP, transitions of the DP class occur in a broad range of systems, including models of catalytic reactions [24], damage spreading transitions [25] and interface growth [26].

At the critical activation rate of the transition, labeled $\lambda_{c}$, the standard CP follows a power-law behavior for long times. When starting from a single active site on an infinite $d$-dimensional lattice, the number of active sites $N$, survival probability $P$, and the the root-mean-square distance of the active sites from the initial active site $R$ scale asymptotically as [18, 19]

[TABLE]

with critical exponents $\eta$, $\delta$, and $1/z$. The standard CP also showcases a phenomenon called duality [27, 28]. When starting from a fully-active lattice with periodic boundary conditions, the density of active sites, $\rho$, is identical to the survival probability of a spreading simulation, started from a single active site, yielding the relationship [18, 19]:

[TABLE]

with $\alpha=\delta$. Duality also implies that the critical infection rate $\lambda_{c}$ is the same for the two initial conditions. Note that the critical exponents above are not independent, as they are related by the hyperscaling relation $\eta=\frac{d}{z}-\delta-\alpha$, where $d$ is the dimension of the underlying lattice [18, 19].

The rest of the paper is organized as follows. The contact process with Hebbian learning is introduced in section I. The impact of Hebbian learning is discussed in section II in 1D and in section III in 2D and above, followed by a Discussion. Detailed analytic and numerical results in 1D are presented in the Appendix A for the imbalance variable, with further results for positive reinforcement in Appendix B.

I Model and Methods

The CP is controlled by an activation rate $\lambda$ and a rate $\mu$ of spontaneous deactivation, which is usually chosen to be one for convenience. In the numerical implementation, we follow the standard simulation prescription by Dickman [29]. An active site $i$ is chosen at random from a list of $N$ active sites. The site transitions to the inactive state with probability $1/(1+\lambda)$ and attempts to activate a neighboring site with probability $\lambda/(1+\lambda)$. In order to approximate a continuous-time Markov chain, each of these steps increases time by $\Delta t=1/N$, so that, on average, each active site acts once per unit of time. In the standard CP, if activation is chosen, a neighboring site is chosen at random with uniform probabilities. If this site is inactive, the activation attempt is deemed successful and the neighbor transitions to the active state. For actual biological, social and physical networks, relationships between sites need not be uniform. To model this, we introduce edge weights, carried in a weight matrix, $W$, where $w_{i,j}$ is the weight of the connection from site $i$ to site $j$. Then, the direction of activation is chosen with probability proportional to the edge weight in that direction. As a simple example, if site $i$ in a 1D chain has chosen to activate a neighbor, it will attempt to activate site $i+1$ with probability $w_{i,i+1}/(w_{i,i-1}+w_{i,i+1})$ and it will attempt to activate site $i-1$ with probability $w_{i,i-1}/(w_{i,i-1}+w_{i,i+1})$.

To model Hebbian learning in the CP, we choose an event-based, multiplicative type of reinforcement rule: the weights are multiplied with a constant factor $r$ after each successful activation event between sites $i$ and $j$. This can be realized according to one of the three possible scenarios: mutual, source or target learning, as illustrated in Fig. 1 and summarized in Table 1. In the mutual learning model, both $w_{i,j}$ and $w_{j,i}$ are multiplied by the reinforcement parameter, $r$. In the source learning model, only the weight $w_{i,j}$ is multiplied by $r$, while in the target learning model only $w_{j,i}$ changes. In all cases, if site $i$ attempts to activate the already-active site $j$, $w_{i,j}$ and $w_{j,i}$ are left unchanged. In line with the usual intuitive picture of Hebbian learning, reinforcement happens in these models only at events when two sites become simultaneously active. However, simultaneous activity alone does not lead to reinforcement, as a direct activation event is needed.

Note that even with the application of a reinforcement rule, the total rate of activations sourcing from a given site is invariant ($\lambda$), as changes in the weights only affect the preferences according to which a site spends it’s contact time. In other words, we have an implicit normalization for the total activation rate and learning induces directional selectivity towards or against some partners. Without some sort of normalization, Hebbian learning is in general unstable as either all activation rates grow until a maximum allowed value or all rates decay to zero [30].

Our implementation of local learning is in stark contrast to past attempts in the literature to capture adaptation in spreading processes. In the traditional framework of adaptive networks [31, 32], the topology of the connections is changing, while in the more recent approach of “awareness” [33, 34], the infection rate changes according to the local or global activity. Hence, our implementation of Hebbian learning is a qualitatively different adaptation mechanism that might in practice appear simultaneously with other modes of adaptation, such as rewiring the topology or adjusting the local infection rates based on awareness.

To investigate the effects of our Hebbian learning rules on the CP, we consider two initial conditions—spreading from a single active site or decaying from a fully-active lattice. For spreading simulations, we monitor the number of active sites $N$, the survival probability $P$, and the mean square radius $R$; while for a fully active initial condition we monitor the density of active sites $\rho$. Initially, $w_{i,j}=1$ for all neighbors $(i,j)$. Each edge weight can then be written as $w_{i,j}=r^{u_{i,j}}$, where $u_{i,j}$ is the number of times the connection $w_{i,j}$ has been reinforced. In the 1D chain case, the quotient $w_{i,i-1}/(w_{i,i-1}+w_{i,i+1})$ can be written as $1/(1+r^{n_{i}})$ where $n_{i}=u_{i,i+1}-u_{i,i-1}$ is the difference in usage between outgoing edges. Likewise, $w_{i,i+1}/(w_{i,i-1}+w_{i,i+1})$ can be written as $1/(1+r^{-n_{i}})$. Hence, in the 1D chain, all relevant information concerning the actions of site $i$ can be stored in the quantity $n_{i}$, which we refer to as the imbalance of the site. With this definition, the imbalance of a site is zero if it has activated both of its neighbors the same number of times, indicating no bias in the preferences for future infection events.

We note that the contact process endowed with the above event-based learning rules is still a Markov process, however, on an extended state space: besides the usual occupation variables of the standard CP, the configurations of the system are also characterized by the additional set $\{u_{i,j}\}$ (or $\{n_{i}\}$ in 1D) which encode the actual weights. This model differs from the standard CP in that it has an infinite number of absorbing configurations differing only in the weight variables. As a consequence, our model is exempted from the scope of the DP conjecture mentioned in the Introduction.

Note that to check the emergent behavior for values of $r<1$, it is also useful to investigate the limiting case $r\to 0$. When tracking the number of times each edge has been used, as $r$ approaches 0, the site will choose to spread activity along the least used edge more often. In the $r\to 0$ limit, activity always spreads along the least-used edges. If there are more than one of such edges, the direction of activity spread is chosen uniformly randomly, like in the standard CP. Similarly, for $r>1$ it is useful to consider the deterministic case with $r\to\infty$, where activity spreads along the most-used edges.

For investigating spreading from a single active site, we chose our lattice size such that the spread never touches the boundary. For example, to run simulations to $t=10^{6}$ for mutual learning in 1D, it was sufficient to have lattice sizes of $L=10^{6}$ for $r\to 0$. For the same model in 2D, $L=2.5\,10^{4}$ was sufficient up to time $t=10^{5}$. We typically collected over $10^{4}$ samples for spreading simulations and around $10^{2}$ for fully infected simulations. In 1D, fully infected simulations had a typical linear size of $L=10^{6}$, and in 2D $L=5\,10^{3}$, but at least $L=2\,10^{3}$ to reduce finite-size effects.

II Negative Reinforcement in one dimension

In one dimension, negative reinforcement ($r<1$) results in qualitatively different behavior in the three different learning modes as illustrated in Fig. 2. Here, we show the density profiles at various times starting from a single active site at the center. In the source learning model, the activity peaks at the center and spreads according to the DP dynamical scaling, as illustrated at the critical point in Table 2. We find that even the scaling function of density profiles—plotted in the inset of Fig. 2—has the same non-Gaussian tail, $\tilde{\rho}(x)\sim e^{-Cx^{\Phi}}$ with $\Phi\approx 2.4$ as observed for the standard CP [20]. The critical mutual learning model has a qualitatively different behavior: the activity peaks at two ballistically advancing front regions. In stark contrast, there is no active bulk phase in the target learning model, with only two narrow fronts propagating outwards ballistically that have a finite lifetime.

In the following, we present detailed scaling results for the critical source and mutual learning models, supported by analytic arguments and further measurements in the Appendix A on the behavior of the imbalance variable. For the target learning model, there is no stationary state with a positive density of active sites, no matter how large $\lambda$ is, as discussed in the Appendix B.

II.1 Source learning model

We have explored several values of the learning rate $r$, leading to qualitatively similar results for all $r<1$ values. Here, as an illustration, we focus on presenting the results at $r=0.8$, shown in Fig. 3, with additional results summarized in Table 2. We find that the critical point of the standard CP is inactive with negative reinforcement, and the critical point is monotonically shifted towards larger values of $\lambda$ as $r$ decreases. In terms of standard observables (Fig. 3), all $r<1$ cases are similar in the sense that within the error they cannot be distinguished from the DP universality in terms of $\alpha$, $z$, and $\eta+\delta$. Following directly from the definitions, $\eta+\delta$ governs the total number of active sites in the surviving runs as $N_{s}\sim t^{\eta+\delta}$. This means that the observables conditioned on survival, $N_{s}$ and $R(t)$ show DP critical behavior, but the survival probability itself is affected by reinforcement, see the altered values of $\delta$ in Table 2. This entails the change of $\eta$, characterizing the number of active sites in all trials. The latter two exponents seem to change hardly with the reinforcement parameter. Such a scenario is not unprecedented for the CP: the surface critical behavior is also known to be characterized by a modified order parameter exponent, while the exponents related to observables conditioned on survival are unaffected [20]. Note that duality is violated in this model, as evident from the different values of $\alpha$ and $\delta$, although the transition happens at the same $\lambda_{c}$ for both initial conditions. We refer to this type of critical behavior with certain exponents governed by DP, while others (here $\delta$ and $\eta$) deviating from DP, as a manifestation of the learning DP (LDP) behavior.

A deeper understanding of this LDP behavior can be obtained through the imbalance variable. As introduced earlier, the imbalance of a site $n_{i}(t)$ is the signed sum of local successful infection events from site $i$, $-1$ meaning left infection, $+1$ meaning right infection. As shown in Fig. 10b, the imbalance distributions tend to a limit distribution at late times, which has Gaussian tails $p_{n}\sim e^{-Cn^{2}}$ with and $r$-dependent constant $C$ for $|n|\gg 1$. This is found to hold irrespective of whether the model is in the active phase or at its critical point. The Gaussian form is in line with our analytic expectations based on a simple mean-field approximation presented in the Appendix A. We note that, the local variable $n_{i}$ at a fixed site has the same limit distribution as the ensemble of local variables, and this limit distribution is symmetric under $n\to-n$. The Gaussian tail of the limit distribution of imbalance variables translates to a log-normal tail of the distribution of instantaneous weights.

II.2 Mutual learning model

Just like in the source learning case, we obtained qualitatively similar results for all values of $r<1$, illustrated here through the case of $r=0.9$, see Fig. 4, with all results summarized in Table 2. As opposed to source learning, the critical activation rate weakly decreases with decreasing $r$. Thus the $r<1$ mutual learning model is active at the critical point of the standard CP, indicating that mutual negative reinforcement promotes infection spreading. The other remarkable difference compared to source learning is that the survival probability $P(t)$ tends to a positive limit at late times even at the critical point. This means that regarding, as usual, the limiting value of $P(t)$ as an order parameter, the transition is discontinuous, with a formally vanishing exponent, $\delta=0$. The other measured critical exponents also differ from those of the source learning model. The dynamical exponent is close to $1$ for all $r<1$, indicating a ballistic spreading of the active zone, in agreement with the density profiles in Fig. 2b. The asymptotic velocity of spreading is found to increase with decreasing $r$. The number of active sites $N$ shows a power-law like scaling with a much larger $\eta$ exponent than for the standard CP or the source learning model. The global density $\rho(t)$, when the process starts from a fully active initial state, decays algebraically in time, and the exponent $\alpha$ is close to the DP value for not too high values of $r$, otherwise it has a somewhat increased value compared to DP. Duality of the DP class is clearly broken—even the corresponding critical exponents are different $\alpha\neq\delta$. Mutual learning in one dimension thus serves as another example of an LDP, a mixture of DP-like and non-DP behavior, although it is clearly different from the behavior of the source learning model.

The inspection of the density profiles shown in Fig. 2b, as well as the measured values of dynamical critical exponents lead us to the following picture of activity spreading at the critical rate. At a given time, the system can be divided in three parts. An intact region, which has not yet been visited by activity, a front region with a high $[O(1)]$ local density of activity, and an interior part in which sites were reactivated many times so that the imbalance distribution is close to its limiting distribution. It is easy to see that the mutual learning rule with $r<1$ facilitates activation processes from the actual front particle toward the intact region. When starting from a single active site, as activity spreads outwards, the edge connecting the newly active site back towards the center of the cluster is weakened. When it is this newly activated site’s turn to activate others, it is significantly more likely to activate outwards, towards inactive sites. This reduces the amount of times activation is attempted but failed, promoting spreading.

This mechanism results in a ballistic drift of the fronts, at which the process is locally supercritical. In the interior part, the process is critical and the local density decays algebraically in time. This ballistic expansion of the critical zone (c.f. the sublinear creeping of the front $R(t)\sim t^{1/z}$ in the standard CP) ensures that the process remains surviving with a finite probability even at the critical point. Specifically, at a fixed site $x$, the local density is $O(1)$ at time $t_{x}\simeq x/v$, i.e. the time at which the front sweeps through site $x$. Here, $v$ denotes the mean asymptotic velocity of the front. After this event, the local density decreases in time asymptotically as $\rho_{x}(t)\sim(t-x/v)^{-\alpha}$. Summing up these local densities in the interior part $[-vt,vt]$, we obtain an algebraic increase of the number of active sites $N(t)\sim t^{1-\alpha}$. Assuming that the bulk criticality belongs to the DP universality class, we obtain $1-\alpha=0.8405$. Numerical estimates presented in Table 2 are somewhat lower; this deviation may be attributed to corrections stemming from the front region which is found to spread out proportionally to $t^{1/z}$ with the DP dynamical exponent $z$ on top of the ballistic advancement (not shown).

As shown in Fig. 11, the distribution of the imbalance values tends to a stationary distribution with an exponential tail, $p_{n}\sim e^{-an}$, with $a=0.87(2)$. This scaling form is in agreement with our mean-field calculation, shown in the Appendix A. Note that the exponential tail of imbalance variable corresponds to a power-law tail of the instantaneous weights.

III Higher dimensional systems

In higher dimensions, for simplicity, we restrict our presentation to the results obtained in the limits $r\to 0$ and $r\to\infty$. These limits are computationally more feasible to study, lead to less severe scaling corrections and are usually expected to be representative of the generic cases with negative ($r<1$) and positive ($r>1$) reinforcement, respectively. As in 1D, the target learning model is an exception at $r\to 0$, as in this limit the heterogeneities become frozen, leading to an emergent quenched disorder, as discussed later in this section.

III.1 Negative reinforcement

In the source learning case, our simulation results are illustrated in Fig. 5 for $r\to 0$, with the estimates of the critical exponents listed in Table 3. For a fully-active initial condition, $\lambda_{c}\approx 1.695$ and $\alpha$ is around the DP value, while spreading simulations indicate a lower critical rate around $\lambda_{2}\approx 1.6768$. In the mutual learning model, we again observe two critical points, as illustrated in Fig. 6, with the exponents presented in Table 3. Our results indicate that the system is moved to a supercritical state when run at the clean critical point but with $r<1$. Both the spreading critical point ($\lambda_{1}$) and the bulk critical point ($\lambda_{2}$) is shifted down for $r<1$ negative reinforcement, just like in the 1D case.

A similar scenario of two phase transitions is observed for the target learning model, as shown in Fig. 7. Here, $\lambda_{1}$ is below while $\lambda_{2}$ is above the clean critical point. For this model with $r\to 0$, a special spontaneous ”freezing” of random directions takes place which is analogous to the extraordinary behavior of the same model in one dimension. Here, the set of least-used (smallest-weight) directions at each site constitutes the “directional configuration” of the system. Occupation variables then follow a simple dynamics that is completely determined by the directional configuration of the system. It is easy to see that there is an infinite number of directional configurations which remain invariant under the dynamics of the model. The only restriction for creating such invariant configurations is non-reciprocity, i.e. whenever site $i$ has a least-used direction to site $j$ then site $j$ must not have a least-used direction to site $i$. A reciprocal situation is unstable against the dynamics and can only happen if the local activity ceases to exist before fixing the direction of the link. Indicating the least-used directions by arrows on links, invariant configurations are constructed by placing arrows (in either direction) on links or not placing any, under the condition that each site has at least one outgoing arrow. According to our numerical observations, the $r\to 0$ target learning process—when started from an isotropic directional configuration—evolves under the dynamics to an invariant directional configuration, which is chosen randomly by the fluctuating dynamics. Obviously, the evolution of directional configuration may also freeze locally in some domain before settling into an invariant configuration by the vanishing of activity in that domain. This quenching phenomenon is similar to what happens in one dimension, although the directional configuration there is trivial due to the 1D topology.

In all cases, our observations imply the existence of two phase transitions in these models, dividing the emergent behavior into three distinct phases of i) no activity below $\lambda_{1}$, ii) a propagating front with no bulk activity for $\lambda_{1}<\lambda<\lambda_{2}$, and iii) an active bulk phase for $\lambda>\lambda_{2}$. As a further evidence for this picture, at $\lambda_{1}$ we eventually observe a fast decay of the local density at the origin as illustrated in Fig. 8, indicating that the bulk is asymptotically inactive at $\lambda_{1}$. The $\lambda_{2}$ location of the bulk critical point is different in all three models, as identified by the critical point of fully infected simulations based on the scaling of the total density $\rho(t)$.

The source and mutual learning models are close to the DP scaling at the bulk critical point, $\lambda_{2}$. In the case of the target learning model, for which spontaneous freezing of directions occurs (cf. the fluctuating weights in the other two models), the exponent $\alpha$ differs significantly from the DP value, although strong corrections may also not be excluded. At the spreading critical point $\lambda_{1}$, we see a marked deviation from the DP expectations, close to the values of the dynamical percolation (DyP) universality class. In this case, target learning is closest to the DyP exponents, but the other two model variants might also show the same behavior (as suggested by $z$ and $\delta$), albeit with stronger corrections (mostly in the case of $\eta$). The DyP is known to capture the phase transition in the susceptible-infected-recovered (SIR) model. Note that the DyP exponents satisfy a modified hyperscaling relation [20], given by $\eta=d/z-2\delta-1$. This hyperscaling relation is also satisfied by our exponents for each of the three learning models at $\lambda_{1}$. Next, we argue that indeed, at $\lambda_{1}$, asymptotically all three learning models are expected to follow a SIR dynamics.

As in the 1D model with mutual learning, we can distinguish three different domains of the spreading system above $\lambda_{1}$. These are the i) ballistically expanding annular front region, where the activity is high, ii) the unvisited domain outside of the front, and iii) the interior part which has already been swept through by the front. Spreading of the front region toward the outer empty region where there is no directional bias is facilitated compared to spreading in the inner domain, in which there is a random directional field for transmitting activity. As a consequence, for a not too high value of $\lambda$, the process is still locally subcritical in the interior, although the front region is steadily expanding outwards. Note that once activity visited a local region, the homogeneity of weights there is irreversibly broken. This suggests an analogy with the dynamics of the SIR model, which differs from the standard CP in that active sites do not return to the inactive state but get removed irrevocably from the system. The analogy is salient: the unvisited (homogeneous) local regions correspond to the susceptible state of the SIR model, whereas the visited (inhomogeneous) local regions correspond to the recovered state. This correspondence is obviously not perfect: in the SIR model, recovered sites can never be reactivated, while in our model this rule is less strict, since activity is able to penetrate into an inhomogeneous region. Nevertheless, as this latter region is subcritical, the penetration length is finite. Hence, concerning large-scale properties at the critical point, this circumstance is irrelevant. Based on this observation, we expect the lower (spreading, $\lambda_{1}$) critical points of our two-dimensional learning models to belong to the DyP universality class of the SIR model.

Thus, the lower phase transition of the model ($\lambda_{1}$) is governed by the criticality of the front region, while the interior of the system is subcritical at this point. It is clear that such a situation can be realized only in dimensions higher than $1$, in which the boundary between the intact and visited domain has a dimension $d-1>0$. On the contrary, in one dimension the criticality of the system cannot be sustained solely by a (zero-dimensional) front. As a consequence, in 1D models, the entire inner region must be critical at the spreading transition. This is in agreement with our numerical observations that our models in one dimension display only a single phase transition.

III.2 Positive reinforcement

As found in Appendix B, in the one-dimensional case within a mean-field treatment, there is generally no active phase with positive reinforcement. This is also the case in higher dimensions. While the mutual and target learning models lead to a quick localization of the activity for all $r>1$, the source learning case can show a qualitatively different behavior in higher dimensions. This generic phenomenon is easier to understand for large $r\to\infty$. In this case, it is easy to see that the direction of the first successful activation from a given site ”freezes” i.e. any further activation from that site can occur only in that direction. This means that learning leads to an emergent, correlated quenched disorder in the system. As a consequence, in spreading simulations, a train of particles is forming, the front particle performing a random walk until it revisits a site that was activated before. In any dimensions, this happens with probability one—after a finite transient—and then the activity can only propagate around the one-dimensional loop. This is a manifestation of the “ant mill” phenomenon [41]. An ant mill emerges when certain ants reinforce their trajectories by pheromones following each other until they die of exhaustion when accidentally forming a loop. The simplest mathematical model considers multiplicative reinforcement in the form of a directed edge reinforced random walk, with rigorous results in all $d\geq 2$ dimensions [42].

Large loops happen with an exponentially small probability, but they are active for an exponential amount of time above the critical point of the one-dimensional totally asymmetric CP with $\lambda^{*}=3.306(2)$ [40]. Hence, above $\lambda^{*}$ we observe a non-universal power-law decay of the survival probability as shown in Fig. 9 in the large $r$ limit, a hallmark of a Griffiths phase [43, 44]. The emergence of a Griffiths phase is expected to be a generic result that holds for all $r>1$ reinforcement rates in all $d\geq 2$ dimensions above $\lambda^{*}$. Note that in one dimension, loops have only two sites, so a Griffiths phase with the same mechanism cannot occur.

Discussion

In this study, we introduced generalizations of the contact process with a local learning rule inspired by Hebbian learning, leading to rich emergent phenomena. We found that local incentives can lead to the opposite global effects, like locally avoiding previously activating or activated neighbors can turn the inactive phase into a globally active phase, and vice versa. For example, positive reinforcement in the mutual learning model is akin to quarantining. In this case, if infection is passed between two sites, they preferentially interact with each other, instead of their other neighbors, leading to activity localization at a finite scale that depends on the value of the reinforcement parameter, $r>1$. such localization hinders activity spreading at the global scale, forcing the system to be in the inactive phase.

In stark contrast to standard critical phenomena, microscopic differences in the learning models can lead to dramatic changes in the emergent behavior, from completely destroying a phase transition to introducing an additional phase transition. For example, in two dimensions, with negative reinforcement, we see different phase transition depending on the initial conditions. At the fully infected (bulk) critical point, we observe a learning DP (LDP) behavior, where some observables follow the DP scaling, while others show a distinct behavior affected by learning. At the spreading critical point, in higher dimensions, our results suggest a critical point in the universality class of dynamical percolation (DyP), as in the SIR model. In this case, learning with negative reinforcement turns the SIS model (without immunization) into an effective SIR model, even though the local infection and healing rates remain unaffected. In the region between the spreading and the bulk critical point, the system is in a phase showcasing a propagating front with an subcritical bulk. It is an interesting future direction to explore if such a phase could be observed in some real-life applications, either in infection, brain activity or population spreading. Although the contact process is a genuinely non-equilibrium model and learning leads to changes in the model parameters coupled to the dynamics, we were also able to obtain informative analytic results, further supporting our simulation observations.

In the target learning model, brain-inspired positive reinforcement leads to a Griffiths phase in higher dimensions. There is increasing evidence that brain dynamics shows hallmarks of being in a Griffiths phase [45, 46, 47, 48]. However, in standard models, a Griffiths phase is either caused by uncorrelated time-invariant (quenched) disorder [43], or a hierarchical-modular architecture of cortical networks [8, 49]. Our results serve as an independent alternative mechanism, a manifestation of the ant mill phenomenon [42], originating from dynamical trapping in activity loops emerging from correlated heterogeneities induced by learning. In the brain, it is believed that being in a Griffiths phase offers functional advantages, in terms of memory and computational capabilities [8]. Signatures of Griffiths phases might be also present in artifical neural networks [50], an exciting future direction to explore in relation to Hebbian learning.

The target learning model also leads to an emergent quenched system in the $r\to 0$ limit, leading to invariant directional configuration(s) in both 1D (totally asymmetric CP [40]) and 2D. Note that in 2D the system is disordered, as there are infinitely many directional configurations depending on the initial condition and the dynamics. The directional configuration also contains correlations resulting from the generating dynamics. Without learning, uncorrelated quenched disorder in the rates leads to a distinct type of critical behavior characterized by ultra-slow activated scaling [43] below four dimensions [37, 51] and on certain types of complex networks [52], in line with the perturbative argument known as the Harris criterion [53]. It is an interesting open question to explore whether such a directional quenched disorder is a relevant perturbation at the critical point, potentially leading to exotic, strong disorder phenomena.

To the best of our knowledge, our work presented here is the first study of Hebbian learning for activity spreading, with several future directions remaining open. Such extensions could consider alternative forms of Hebbian learning, either incorporating information on the duration when the sites are simultaneously active, or different normalization schemes. For example, in the case of positive reinforcement, brain dynamics would suggest the study of the potentially stabilizing impact of inherently decaying activation rates [30]. As real-life infection spreading is complicated, combinations with other mechanisms, like the simultaneous impact of Hebbian learning and awareness [33, 34] could lead to more realistic models.

Although the CP is equivalent to the SIS model on a lattice, this is no longer true on complex networks. In stark contrast to the CP [54], in the SIS model there is no inactive phase on scale-free networks with degree exponent $\gamma>2$ [55, 56, 57]. Hence, positive reinforcement might not completely eliminate the active phase. Instead, a phase transition might happen between an inactive and an active phase on such networks. Another direction of study is unveiling the condition under which a similar reinforcement-driven Griffiths phase occurs on complex networks through the ant mill phenomenon [42].

Acknowledgements.

This work was supported by the National Research, Development and Innovation Office NKFIH under Grant No. K146736. The work of I.A.K. was supported by the National Science Foundation under Grant No. PHY-2310706 of the QIS program in the Division of Physics, by the Domus Hungary Scholarship of the Hungarian Academy of Sciences, and by the Baker Faculty Grant of the Weinberg College of Arts and Sciences, Northwestern University. W.T.E. was supported by the Northwestern Office of Undergraduate Research’s Summer Undergraduate Research Grant (SURG) and the Baker Program in Undergraduate Research. We thank L. V. Luzzatto, N. Love, and J. Barrera Lopez for useful discussions, and M. M. Driscoll, B. Oborny and G. Ódor for useful comments.

Appendix A Imbalance in one dimension

A.1 Standard contact process

For the standard CP of $r=1$ the rates are constant, i.e. independent from the imbalance variables $n_{i}$. According to our numerical observations, the distribution $p_{n}$ gets broader algebraically as $t^{\kappa}$. In the active phase, the correlation time $\tau$ is finite. Therefore, the imbalance variables $n_{i}(t_{m})$ sampled at times $t_{m}=mT$, where $T\gg\tau$ and $m=0,1,2,\dots$, are expected to have uncorrelated increments, leading to $\kappa=1/2$. This is indeed in agreement with our simulations (not shown). At the critical point, however, we have found that the scaling of the distribution $p_{n}$ with time is characterized by a smaller exponent $\kappa=0.43(1)$, see Fig. 10a. Besides possible temporal correlations in the process $n_{i}(t)$, there is also an obvious reason for the difference between the exponent $\kappa$ in the active phase and at the critical point. Clearly, the change of imbalance variables is triggered by local activity and the ’clock’ of this process is ticking only if the corresponding site is active, otherwise it is still. In the active phase, the mean local activity is constant in time in the steady state. As opposed to this, at the critical point, it decreases as $\rho\sim t^{-\delta}$. As a consequence, the time $t^{\prime}$ relevant for the process of imbalance variables at the critical point is related to the real time $t$ as $t^{\prime}\sim\int^{t}\rho(t^{{}^{\prime\prime}})dt^{{}^{\prime\prime}}\sim t^{1-\delta}$. Ignoring temporal correlations, the width of the distribution would then increase proportionally to $t^{\prime 1/2}\sim t^{(1-\delta)/2}$, yielding

[TABLE]

at the critical point. This is compatible with the numerical estimate, suggesting that temporal correlations play a weak role, if any.

A.2 Source learning model

A.2.1 Dynamics

The imbalance $n$ keeps track of the cumulative number of successful activation events from a given site, which counts the activations to the right (left) with a positive (negative) sign. The activation rate to the right is $r_{n}=\lambda\frac{r^{n}}{1+r^{n}}$, whereas the rate to the left $l_{n}$ is fixed by $l_{n}+r_{n}=\lambda$. Thus, $r_{n}=l_{-n}$. In the source learning model, the imbalance of the source site is increased (decreased) by $1$ after each successful activation to the right (left). We will denote the state of an active (inactive) site with imbalance $n$ by $A_{n}$ ($I_{n}$). The probability of finding a site in state $A_{n}$ and $I_{n}$ is given by $p_{n}$ and $q_{n}$, respectively. The probability that the site is active is $p\equiv\sum_{n}p_{n}$, and we also introduce $q\equiv 1-p$.

Let us consider a mean-field approximation of model, in which the local states of different sites are assumed to be uncorrelated. For a technical convenience, we will restrict ourselves in the subsequent analysis to the domain $r>0$. The terms contributing to the change rate of $p_{n}$ at a marked site are then the following. i) The loss term by spontaneous deactivations: $-\mu p_{n}$. ii) The loss by activation events from the marked site to the right, $-r_{n}p_{n}q$, and to the left, $-l_{n}p_{n}q$. In total $-\lambda p_{n}q$. iii) The gain by activation events from the marked site to the right $r_{n-1}p_{n-1}q$, and to the left, $l_{n+1}p_{n+1}q$. iv) The gain by the activation of the marked site from the left, $(\sum_{m}r_{m}p_{m})q_{n}$, and from the right, $(\sum_{m}l_{m}p_{m})q_{n}$. In total $\lambda pq_{n}$. The rate of change of $p_{n}$ is then

[TABLE]

Summing over $n$, we obtain

[TABLE]

which coincides with the mean-field description of the standard CP. Thus, within this approximation, the reinforcement does not shift the location of the critical point, $\lambda_{c}/\mu=1$.

Now we turn to the rate of change of $q_{n}$. The corresponding terms are the following. i) There is gain by spontaneous deactivations, $\mu p_{n}$. ii) There is a loss term by activations of the marked site from the left, $-(\sum_{m}r_{m}p_{m})q_{n}$ and from the right, $-(\sum_{m}l_{m}p_{m})q_{n}$. In total $-\lambda pq_{n}$. Hence, we have

[TABLE]

A.2.2 Stationary state

In a weak sense, $\dot{p}=0$, the model has a non-trivial ($p>0$) stationary state above the threshold, $\lambda/\mu>1$, for any $r>0$. However, non-trivial stationarity in a stronger sense, $\dot{p}_{n}=\dot{q}_{n}=0$, exists only for $0<r<1$ i.e. for negative reinforcement.

To see this, let us assume that the process is stationary in the weak sense. Then Eq. (7) implies that the probability of finding a state in an active state is constant, $p=1-\mu/\lambda$. In this case, Eqs. (6) and (8) can be regarded as the master equation of a continuous-time random walk on a comb-like lattice. The allowed transitions are nearest-neighbor hops in the set of active states, $A_{n}\to A_{n+1}$ and $A_{n}\to A_{n-1}$, with rates $r_{n}q$ and $l_{n}q$, respectively. As well as hops to “dangling” states $I_{n}$ back and forth, $A_{n}\to I_{n}$ and $I_{n}\to A_{n}$, with rates $\mu$ and $\lambda p$, respectively. It is easy to see that this random walk will be recurrent only if $0<r<1$, in which case it has a stationary state in the strong sense, $\dot{p}_{n}=\dot{q}_{n}=0$. This leads to

[TABLE]

and

[TABLE]

As the transition graph is acyclic, detailed balance holds in the stationary state, $r_{n}p_{n}=l_{n+1}p_{n+1}$. The stationary probabilities can then be written in terms of a potential $U_{n}$, as $p_{n}=\mathcal{N}e^{-U_{n}}$, where the potential differences are fixed by $U_{n+1}-U_{n}=\ln\frac{l_{n+1}}{r_{n}}$ and $\mathcal{N}$ denotes a normalization constant. Let us now consider the large-$n$ asymptotics of stationary probabilities $p_{n}$ in the case $0<r<1$. For $n\gg 1$, we have $r_{n}\simeq\lambda r^{n}$ and $l_{n}\simeq\lambda$. The potential difference is then $U_{n+1}-U_{n}\simeq n\ln(1/r)$. The potential for large $n$ is thus $U_{n}\simeq n^{2}\frac{\ln(1/r)}{2}$. The stationary probability distribution of the imbalance has therefore the tail

[TABLE]

This Gaussian tail of the imbalance distribution is also observed in simulations, see Fig. 10. Note that the characteristic width of the imbalance distribution is $n^{*}\sim[\ln(1/r)]^{-1/2}$, which is diverging as $r\to 1$. Nevertheless, the asymmetry between activation rates to the right and to the left, which can be conveniently measured by $a_{n}\equiv\ln(r_{n}/l_{n})$, has a characteristic value $a^{*}\sim\sqrt{\ln(1/r)}$, which is vanishing as $r\to 1$.

In the case of a non-negative reinforcement, $r\geq 1$, the random walk is no longer stable, and there exists no stationary state. At $r=1$, the walker spreads out diffusively, $\langle n^{2}\rangle\sim t$, while for $r>1$, the spreading is ballistic, $\langle n^{2}\rangle\sim t^{2}$. For a weak positive reinforcement, $r\gtrsim 1$, the dynamics of the imbalance variable will display a crossover phenomenon from the non-reinforced behavior at short times, where the random walk is recurrent to the positively reinforced case at late times, where the random walk is transient. The crossover time scale $\tau$ will diverge as $r\to 1+$. The form of the divergence within the mean-field approximation can be obtained as follows. The potential for $n\geq 1$ can be evaluated to be $U_{n}-U_{0}=-\frac{n(n+1)}{2}\ln r-\ln(\frac{1+r^{n}}{2})$. If $r$ is close to $1$, we have $U_{n}-U_{0}=-\frac{n^{2}}{2}(r-1)+O(n^{2}(r-1)^{2})$. The crossover from recurrence to transience occurs when the walk reaches a point of no return $n^{*}$ which is marked by the condition $U_{n^{*}}-U_{0}\sim O(1)$. This scales as $n^{*}\sim(r-1)^{-1/2}$, and, since the random walk is essentially unbiased for $|n|\ll n^{*}$, we find for the crossover time scale

[TABLE]

as $r\to 1+$, in line with our numerical results in Fig. 12.

A.3 Mutual learning model

A.3.1 Dynamics

The mutual learning model differs from the source learning model in that, after a successful activation event, also the imbalance of the target site is changed by $1$, with an opposite sign compared to that at the source site. Concerning the change rate of $p_{n}$, this extension modifies only term iv) in the previous section. Here, the gain by activation events of the marked site from the left will be $(\sum_{m}r_{m}p_{m})q_{n+1}$, while from the right $(\sum_{m}l_{m}p_{m})q_{n-1}$. The rate of change of $p_{n}$ is thus governed by

[TABLE]

where $\overline{rp}\equiv\frac{1}{p}\sum_{m}r_{m}p_{m}$ and $\overline{lp}\equiv\frac{1}{p}\sum_{m}l_{m}p_{m}$. Summing over $n$ leads again to Eq. (7), thus $\lambda_{c}/\mu=1$ holds also in the mean-field approximation of the mutual learning model. The rate of change of $q_{n}$ is the same as in the source learning model, see Eq. (8).

A.3.2 Stationary state

Similarly to the source learning model, there is a non-trivial stationary state above the threshold in the weak sense for any $r>0$. In the stronger sense, a non-trivial stationary state exists again only for $0<r<1$.

Let us again assume weak stationarity, and for the sake of simplicity, let us also assume that the initial state is symmetric under $n\to-n$. Then it will remain so at any time, and, as a consequence $\overline{lp}=\overline{rp}=\lambda/2$, i.e. $\overline{lp}$ and $\overline{rp}$ are constant in time. Eqs. (13) and (8) can then be interpreted as a master equation of a random walk on a ladder with the following transitions. $A_{n}\to A_{n-1}$ and $A_{n}\to A_{n+1}$ with rates $l_{n}q$ and $r_{n}q$, respectively. Furthermore, hopping to the other leg of the ladder, $A_{n}\to I_{n}$, with rate $\mu$, and hopping back diagonally as $I_{n}\to A_{n\pm 1}$ with rate $\frac{\lambda}{2}p$. As in the source learning case, this random walk will be recurrent and have a stationary state only if $0<r<1$. Setting $\dot{p}_{n}=\dot{q}_{n}=0$ in the dynamical equations, we obtain Eq. (9) and

[TABLE]

Thus, the stationary probabilities $p_{n}$ are, up to a normalization constant, the stationary probabilities of a random walk on the integers with hopping rates $\omega_{n,n+1}=r_{n}+\frac{\lambda}{2}$ and $\omega_{n,n-1}=l_{n}+\frac{\lambda}{2}$. Let us consider the asymptotics of the imbalance distribution for $0<r<1$. The potential difference of the corresponding random walk for large $n$ is given by $U_{n+1}-U_{n}=\ln\frac{\omega_{n+1,n}}{\omega_{n,n+1}}\simeq\ln 3$. Thus, $U_{n}\simeq n\ln 3$, and the tail of the stationary imbalance distribution is exponential:

[TABLE]

As discussed in the main text, this functional form is again in agreement with numerical simulations, although with a modified decay rate, see Fig. 11.

In the case of non-negative reinforcement, $r\geq 1$, the random walk process is unstable and $\langle n^{2}\rangle$ increases unboundedly, in the same way as in the source learning model. The divergence of the time scale as $r\to 1+$, can be obtained in the same way as for the source learning model. Here, we one has to consider that although the potential varies linearly for $|n|\gg(r-1)^{-1}$, the point of no return is closer to the origin, scaling as $n^{*}\sim(r-1)^{-1/2}$. Indeed, for small $r-1$, the potential is parabolic around the origin, $U_{n}-U_{0}=-\frac{n^{2}}{4}(r-1)+O((r-1)^{2}n^{3})$. This leads to the same form of the divergence of time scale as found in Eq. (12) for the mutual learning model, see Fig. 12.

A.4 Target learning model

In the target learning model, after a successful activation event, the imbalance of the target site is updated (in the same way as for the mutual learning model) but the imbalance of the source site is left unchanged. Within a mean-field approximation which ignores correlations between occupation and imbalance variables at different sites, the rate of change of $p_{n}$ reads in one dimension as

[TABLE]

while that of $q_{n}$ is the same as given in Eq. (8). Assuming weak stationarity and symmetry under $n\to-n$, these equations describe a random walk on a ladder with vertical transitions $A_{n}\to I_{n}$ with rate $\mu$, and diagonal jumps $I_{n}\to A_{n\pm 1}$ with rate $\lambda p/2$. This leads to an unbounded diffusive broadening of the imbalance variable as $\langle n^{2}\rangle\sim t$, irrespective of the value of the reinforcement parameter $r$. Therefore, there is no stationary state in this case.

Appendix B Positive reinforcement in one dimension

For positive reinforcement with $r>1$ there is no active phase, as the active sites form isolated localized clusters that eventually inactivate without triggering a macroscopic active phase. The initial state $n_{i}=0$ becomes unstable and the local imbalance variables start to drift either in the positive or in the negative direction, which is singled out randomly by stochastic fluctuations. The return to the origin ($n_{i}=0$) becomes more and more unlikely with increasing $|n_{i}|$, thus the symmetry of local variables present for $r<1$ is spontaneously broken here. Evidently, an active stationary state would then lead to extreme random imbalance $|n_{i}|\to\infty$, which, however, cannot sustain an active steady state. As a consequence, a non-trivial steady state with a positive density cannot exist in the case $r>1$, which is in agreement with our numerical observations. The state of the system when started from a fully active state with zero imbalance will therefore tend to the absorbing state at late times. This process occurs on a finite time scale, the characteristic lifetime $\tau$, which depends on $r$ and expected to diverge as $r\to 1+$. According to mean-field theory presented in the Appendix A, the divergence is of the form $\tau\sim 1/(r-1)$ above the critical activation rate of the standard CP ($r=1$). Numerical results presented in Fig. 12 indeed show an abrupt decline in the time dependence of density at a characteristic time diverging in accordance with the prediction of mean-field theory. At the critical activation rate of the standard CP, $\lambda=\lambda_{c}$, this is to be corrected by the replacement $\tau\to\tau^{1-\delta}$, yielding $\tau\sim(r-1)^{-1/(1-\delta)}$, see Sec. A.1. Numerical results shown in Fig. 12 are agreement with this prediction. For mutual learning the only difference is that at $\lambda=\lambda_{c}$, we find a slight deviation from the expectation $\tau\sim(r-1)^{-1/(1-\delta)}$, see Fig. 12. The optimal collapse (not shown) is obtained by an exponent $1.25$ instead of $1/(1-\delta)=1.1897$.