#  Slow relaxation and aging in the model of randomly connected cycles network

S. Reich, S. Maoz, Y. Kaplan, H. Rappeport, N.Q. Balaban, and O. Agam  The
Racah Institute of Physics, The Hebrew University of Jerusalem, 9190401,
Israel

###### Abstract

We propose a statistical model of a large random network with high
connectivity in order to describe the behavior of E. coli cells after exposure
to acute stress. The building blocks of this network are feedback cycles
typical of the genetic and metabolic networks of a cell. Each node on the
cycles is a spin degree of freedom representing a component in the cell’s
network that can be in one of two states - active or inactive. The cycles are
interconnected by regulation or by the exchange of metabolites. Stress is
realized by an external magnetic field that drives the nodes into an inactive
state, and the time the magnetization passes zero value for the first time
represents the first division event of the cell after the stress period. The
numerical and analytical solutions for this first passage problem reproduce
the aging dynamics observed in the experimental data.

###### pacs:

87.17.-d Cell processes
05.40.-a Fluctuation phenomena, random processes, noise, and Brownian motion
75.10.Nr Spin-glass and other random models

##  I Introduction

The behavior of the cellular network of molecular components in living cells,
from a physicist’s viewpoint, is an intriguing problem [1]. On one hand, a
cell can be viewed as a large computer, evolved by evolution over billions of
years, and programmed to deal with various environmental changes. In
particular, the cell’s network typically has built-in adaptation mechanisms
that allows it to withstand moderate starvation conditions by taking the cell
into a new state that preserves its vitality. On the other hand, the cell is
also a complex physical system with intricate interactions and collective
behavior amenable to some sort of statistical description.

Most experimental and theoretical studies of cells exposed to stress
conditions are focused on the regulatory regime. These aim to decipher the
various layers of the cellular network (metabolism, gene regulation, etc.) by
disentangling its various pathways, either by focusing on specific small
modules [2] or in similarity to the deciphering of a computer’s blueprint [3]
. However, this reductionist approach is limited because the cell’s cellular
network is incredibly complicated and strongly intertwined.

Nevertheless, several approaches have been proposed to model the cell’s
network using statistical models [4, 5, 6], and using similar approaches as
the ones used in neural network theory [7]. The observation that exposing
cells to unforeseen challenges (by genetic engineering) [8, 9] accentuates the
statistical nature of their behavior as complex systems [10] suggests an
alternative approach to advance our understanding of these complicated
systems. In a recent work [11], our goal was to study the cell in a regime
beyond its ability to respond to environmental conditions, under perturbations
that drive the cell to the edge of death but do not kill, and where it behaves
as a random complex system. In this regime, the cell’s dynamics are dictated
by the global network architecture, which is essentially random. Hence, it can
be described by a relatively simple statistical model. Such a model can be
used as a starting point for developing an increasingly improved description
of the cellular dynamics. The premise of this approach is that the statistical
behavior of the cell holds valuable information about the functionality and
the architecture of the cellular network.

In particular, recent experimental data [11] show that the recovery time of E.
coli cells, after their growth was arrested by acute stress, exhibits
distinctive statistical features typical of complex physical systems. Acute
stress can be achieved, for instance, by adding serine hydroxamate (SHX) to
exponentially growing cultures. SHX induces an artificial starvation for the
amino acid serine and results in growth arrest [12]. In the experiment, cells
are exposed to SHX for a duration $t_{w}$ (waiting time), after which it is
washed out, and the stress period ends. Next, the single-cell lag time, i.e.
the time between the end of the stress and the first post-stress division, is
recorded, and the recovery time distribution is measured from a large ensemble
of cells.

Figure 1:  The experimental result of 1$-$cdf associated with the distribution
of first division time of cells after various starvation periods $t_{w}$.
(Adapted from [11])

The experimental data is conveniently presented by plotting 1$-$cdf (cdf
stands for cumulative distribution function), i.e., the fraction of bacteria
still in the lag phase. A log-log plot of 1$-$cdf as a function of time, for
various values of the waiting time, $t_{w}$, is presented in Fig. 1 [11]. This
figure highlights three distinctive features of the system dynamics: First,
the lag-time distribution depends on the waiting time; namely, it features a
memory effect. Second, the tail of the distribution function exhibits a slow
(approximately) power-law decay. Third, the distribution function saturates at
long waiting times. The same features are common to finite physical systems
that exhibit aging [13], such as amorphous polymers [14, 15], stretched DNA
[16], paper crumpling [17, 18], colloidal solutions [19, 20], spin glasses
[21, 22, 23, 24], supercooled liquids [25], and coulomb glasses [26].

We reiterate that aging behavior only appears when the cell experiences acute
stress. It is not manifested when the stress is weak or gradual (e.g., when
starvation is by natural depletion of nutrients or by gradual application of
SHX). In this case, the cells adapt by upregulating intracellular production,
increasing import, or switching to a different metabolism that prepares cells
for survival and regrowth when nutrients become available. With gradual
stress, 1$-$cdf is practically independent of $t_{w}$ and decays exponentially
in time [11].

In this work we construct and solve a minimal toy model that reproduces,
qualitatively, the results shown in Fig. 1. The general approach is similar to
that of Random Matrix Theory for chaotic quantum systems [27], where a system
with complex interactions and chaotic dynamics, such as atomic nuclei, is
regarded as a ’black box’ with a Hamiltonian matrix drawn from a Gaussian
random distribution. This approach reproduces correctly the universal
statistical features of the energy spectra and eigenfunctions of quantum
chaotic systems in the proper energy regime.

In the same spirit, we consider the cell to be a ‘black box’ with unknown and
complicated interactions among its components (proteins, metabolites, enzymes,
RNA, etc.). To simplify the problem, we assume the following: (a) Each cell
component is described by a boolean variable, i.e., a spin, accounting for
only two possible states of the component: active or inactive. (b) The time
evolution is discrete, and the state of the spins at a given time is
determined, deterministically, by the state of the system at its previous time
step (we use zero-temperature Glauber dynamics [28]). (c) The web of
interconnections among the spins (i.e., the cell components) is entirely
random.

The most natural candidate for a description of this type is the Sherrington
Kickpartick model [29] with zero-temperature Glauber dynamics. This remarkably
simple archetype model of glasses exhibits aging, and slow relaxation [30, 31,
32, 33, 34, 35, 36, 37, 38, 39, 40] as we require. However, being Hamiltonian,
the SK model features a reciprocity property not shared by the cellular
network - the coupling constants between the $i$-th spin and the $j$-th spin
are symmetric, $J_{ij}=J_{ji}$. However, living cells lack this symmetry
because, for example, enzymes catalyze the synthesis of products but not vice
versa.

The asymmetric Sherrington Kickpartick model, in which all coupling constants,
$J_{ij}$, are statistically independent, is also ruled out because it
exhibits, essentially, an instantaneous decay of correlations [41, 42]. Thus,
some elements of the cell’s network architecture should be retained to obtain
slow relaxation and aging from an asymmetric spin network.

Here we propose a model that highlights the metabolic cycles structure of the
cell by choosing its primary constituents to be closed cycles of spins rather
than single spins. The sizes of these loops span over a large range of
timescales, from very short timescales to long ones such as those required for
protein degradation (2 min to 50 hours in growth-arrested cells [43]). There
are indications that these time delays are power-law distributed.

Thus the toy model proposed here consists of a collection of closed spin
chains (representing, for example, different metabolic cycles, or closed
regulatory feedbacks) randomly interconnected via a single spin on each cycle.
The dynamics within each spin chain is a simple shift, and the cycle’s length,
$L$ is a random variable chosen from a power-law distribution with power
$-\alpha$, i.e.

$P(L)=\nu/L^{\alpha};~{}~{}~{}~{}L_{\min}\leq L\leq L_{\max}.$ (1)

Here $L_{\min}$ and $L_{\max}$ are the minimal and maximal endpoints of the
interval on which the distribution is defined, while $\nu$ is the
normalization constant. This model, which we call the Randomly Connected
Cycles Network (RCCN), is illustrated in Fig. 2

In this model, the application of external stress is realized by a magnetic
field that polarizes the spins and sets them into an “OFF” position that
represents an inactive state of the corresponding cell’s component. Soon after
the magnetic field is turned off, the magnetization of the system decays back
to its normal state. The first time it crosses zero is assumed to reflect the
time point when the cell returns to its normal state. It is associated with
the first division time of the cell after starvation, which is the quantity
that is measured experimentally. Thus, the primary task of this work is to
calculate the first passage distribution of the magnetization and compare it
to the first division time distribution of the cells after the starvation
period.

The intention of this paper is as follows: First, we write the equations that
describe the dynamics of the RCCN model. Then we present the results of its
numerical simulations that reproduced, qualitatively, the experimental data
shown in Fig. 1. Next, we analyze the system within a mean-field framework and
use the results of this analysis to develop a simple phenomenological
description of the problem. The latter will allow us to characterize the
distribution function of the total magnetization of the system and obtain an
approximate analytical solution for the first passage problem.

Figure 2: An illustration of the Randomly Connected Cycles Network (RCCN)
model. Each site in the network represents a cellular component, and any red
arrow may assume two positions (up and down), denoting the corresponding
component state - active or inactive. The black arrows stand for a strong
coupling that forms the cycles, while gray arrows represent weak coupling
between different cycles. The coupling is via a single site in each cycle.

##  II Definitions and basic equations

Consider a collection of $N$ closed spin chains that form cycles denoted by
the index $i$. Each cycle contains a random number of spins, $L_{i}$, which we
call the cycle length, and each spin may assume one of two values, $\pm 1$ .
We denote by $s_{i}^{(k)}$ the $k$-th spin on the $i$-th cycle and choose
$s_{i}^{(0)}$ to be the spins that form couplings among all cycles (one spin
in each cycle).

Within each cycle the time evolution is a simple shift dynamics,

$s_{i}^{(k)}(t+1)=s_{i}^{(k-1)}(t),~{}~{}k=1,2,\cdots,L_{i}-1$ (2a) while the
connecting spins satisfy the equation
$s_{i}^{(0)}(t+1)=\\!\mbox{sign}\left[s_{i}^{(L_{i}\\!-\\!1)}(t)\\!+\\!\sum_{j\neq
i}J_{ij}s_{j}^{(0)}(t)\\!+\\!h(t)\right].$ (2b)

Here $h(t)$ is the external magnetic field, and $J_{ij}$ are random coupling
constants with zero mean and a constant variance inversely proportional to the
number of cycles, thus

$\langle J_{ij}\rangle=0,~{}~{}\mbox{and}~{}~{}~{}\langle
J_{ij}^{2}\rangle=\frac{\gamma^{2}}{N},$ (3)

where $\gamma$ is a constant that characterizes the coupling strength. The
magnetization of a cycle is defined to be:

$M_{i}(t)=\frac{1}{L_{i}}\sum_{k=0}^{L_{i}-1}s_{i}^{(k)}(t),$ (4)

and the total magnetization of the system is

$M(t)=\frac{1}{N}\sum_{i=1}^{N}M_{i}(t).$ (5)

In a large system with sufficiently strong coupling among the cycles,
($\gamma\gtrsim 1.3$), for a typical realization of $J_{ij}$, and random
initial conditions of the spins, the magnetization, $M(t),$ behaves like a
random process. The magnetic field application increases its average value,
which decays back to zero after the magnetic field is turned off.

The main quantity of interest is the distribution of times, $t$, of first
crossing, $M(t)=0$, after magnetic field application. This distribution,
denoted by $\rho(t)$, is assumed to reflect the distribution of the first cell
division time after the starvation period. In order to avoid data binning, it
is convenient to present the results using the survival probability,

$S(t)=\int_{t}^{\infty}dt^{\prime}\rho(t^{\prime}),$ (6)

also denoted by 1$-$cdf.

Figure 3:  The magnetization survival probability obtained from simulations of
the RCCN model. The black dashed lines are the approximate solution given by
formula (46). The dependence of these graphs on the starvation period $t_{w}$
reflects aging behavior, which saturates for large values of $t_{w}$.

##  III Numerical simulations.

The simulation procedure employed for the study of the RCCN runs along with
the following steps: First, we choose the lengths of the cycles from the
distribution (1). Then we set the coupling constants $J_{ij}$ to be
independent random variables from a normal distribution with mean and variance
given by Eq. (3). Next, Eqs. (2) are iterated for several (2000) time steps in
order to relax the system to some typical state. A constant magnetic field is
applied for a period $t_{w}$, after which we continue to iterate Eqs. (2) and
record the first time $M(t)$ crosses zero. The statistics are gathered from
different simulations in a similar amount to the number of cells monitored in
the experiment, and the first crossing time data is used to construct the
survival probability.

The choice of different $J_{ij}$ at each run implies that the survival
probability results from ensemble averaging rather than averaging over initial
conditions. Because the system is finite (about 400 cycles), by this way, we
avoid the problem of choosing some non-typical realization of the coupling
constants.

In Fig. 3, we present the results of numerical simulations for the survival
probability for various stress periods, $t_{w}$. These are obtained for a
system with $2^{14}$ spins, $\gamma=3/2$, $\alpha=3/2$, $L_{\max}=2500$,
$L_{\min}=1$, and $h=0.8$ during the stress period. These results are in
qualitative agreement with the experimental data shown in Fig. 1 as they
exhibit the following three features: (a) An approximate power-law decay
within a wide time interval. (b) The delay in the relaxation increases with
the stress period, $t_{w}$. (c) The behavior of the system exhibits saturation
at long stress periods.

Figure 4: A single cycle of the system. In the mean field approach, the input
from all other cycles is assumed to be noise whose autocorrelation function is
determined self consistently.

##  IV Mean Field Approach

The mean-field approach to the RCCN model is based on the assumption that in a
large enough system where all cycles are interconnected, the input signal to a
given cycle is a noise term whose statistical properties can be determined
self consistently.

Consider the $i$-th cycle of the system as illustrated in Fig. 4. Here the
input signal from all the other cycles is

$\zeta_{i}(t)=\sum_{j\neq i}J_{ij}s_{j}^{(0)}(t).$ (7)

When the system is large enough, correlations between cycles are small (see
Appendix A) and the input signal is approximately a Gaussian noise with zero
mean. Correlations can be determined self consistently from the spins
autocorrelation function:

$\langle\zeta_{i}(t)\zeta_{i}(t^{\prime})\rangle=\gamma^{2}\left\langle
s^{(0)}_{k}(t)s_{k}^{(0)}(t^{\prime})\right\rangle,$ (8)

where averaging is over the ensemble $J_{ij}$. In a large enough system,
ensemble averaging is the same as averaging over the system’s initial
conditions. This ergodicity property is discussed in Appendix A.

From here on, to shorten the notations we shall suppress the cycle index, and
denote time using a subscript. With these notations Eqs. (2) for the
connecting spin in a cycle of length $L$ reduce to

$s_{t+1}=\mbox{sign}(s_{t-L+1}+\zeta_{t}+h_{t}).$ (9)

It is instructive to first consider the case where all cycles are of the same
length, $L=1$, and zero magnetic field. Let $p(t)$ be the probability that
$s_{t}=1,$ and assume that the initial condition, $p(0)$, is known. The random
noise drives the system into the fixed point $p(t)\to 1/2$. Yet, there will be
no change in the spin state until time $t$ if the noise is within the range
$-1<\zeta_{\nu}<1$ for all time steps $\nu=1,2,\cdots t-1$. Therefore the
probability $p(t)$ evolves according to the following equation:

$p(t)-\frac{1}{2}=V_{t}\left[p(0)-\frac{1}{2}\right],$ (10)

where, for $t\geq 1$,

$V_{t}=\idotsint\limits_{-1<\zeta_{\nu}<1}~{}~{}\prod_{\nu=0}^{t-1}d\zeta_{\nu}~{}f_{t}(\zeta_{0},\zeta_{1}\cdots\zeta_{t-1})$
(11)

is the probability for no spin-flip for any initial state of the spin within
the time interval $(0,t)$. Here $f_{t}(\zeta_{0},\zeta_{1}\cdots\zeta_{t-1})$
is the joint probability distribution function of the noise at $t$ consecutive
time steps. A formal derivation of this result can be found in Appendix B.

Let us assume that the input noise features an exponential decay of
correlation:

$\langle\zeta_{t}\zeta_{t^{\prime}}\rangle=\gamma^{2}\epsilon^{|t-t^{\prime}|},$
(12)

where $|\epsilon|<1$ is, at the moment, an unknown parameter that determines
the time scale of the decay of correlations.

Inverting the covariance matrix of the noise (12) results in a tridiagonal
matrix that allows one to write the joint distribution function of the noise
in the form

$f_{t}(\zeta_{0},\zeta_{1}\cdots\zeta_{t-1})=Z_{t}\exp\left(-\frac{1}{2}\Phi_{t}\right),$
(13a) where $Z_{t}$ is the normalization factor, and
$\Phi_{t}=\frac{1}{\gamma^{2}(1-\epsilon^{2})}\sum_{\nu=0}^{t-2}\left(\zeta_{\nu}-\epsilon\zeta_{\nu+1}\right)^{2}+\frac{\zeta_{t-1}^{2}}{\gamma^{2}}.$
(13b)

The integral (11) can be evaluated in the limit $\gamma\gg 1$, $|\epsilon|\ll
1$ and $1\leq t\ll\gamma^{2}/\epsilon^{2}$ giving

$V_{t}\simeq\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)\mbox{erf}^{t-1}\left(\frac{1}{\sqrt{2}\gamma\sqrt{1-\epsilon^{2}}}\right),$
(14)

where erf$(x)$ is the error function.

To identify the self-consistent equation for $\epsilon,$ notice that $p(t)$
can be interpreted as the conditional probability of $s_{t}$ given $s_{0}$.
Thus, from $\langle s_{t}\rangle=2p(t)-1$ it follows that $\langle
s_{t}s_{0}\rangle=V_{t}$. Substituting this correlation in Eq. (8) and using
Eq.(12) we obtain the self consistent equation for $\epsilon$:

$\epsilon\simeq-\frac{1}{t}\log
V_{t}=-\log\left[\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma\sqrt{1-\epsilon^{2}}}\right)\right]$
(15)

which for large $\gamma$ gives $\epsilon\simeq(2/\pi)^{1/2}/\gamma$. Thus the
decay time of correlations is

$\tau_{1}=-1/\log\epsilon.$ (16)

In the case where all cycles are of length $L>1,$ the dynamics are obtained by
rescaling of time, and the relaxation time in this case is

$\tau_{L}=L\tau_{1}.$ (17)

However, the situation becomes more complicated when the lengths of the cycles
are distributed according to the power-law given by Eq. (1). To reveal the
behavior of the noise, in this case, one can try to solve the problem by
iteration, starting from an initial approximation where the autocorrelation
function of a spin in a cycle of length $L$ is $\exp(-t/\tau_{L})$.
Substituting this approximation in the right hand side of Eq. (8) and taking
the average over the cycles length distribution yields a power-law decay of
the noise correlations $\langle\zeta_{0}\zeta_{t}\rangle\sim t^{1-\alpha}$
within the range $L_{\min}\ll t/\tau_{1}\ll L_{\max}$. This correlation can be
used to construct the first approximation to the noise distribution function,
$f_{t}(\zeta_{0},\zeta_{1}\cdots\zeta_{t-1})$. With this approximation, the
integral (11) can be calculated in order to obtain the next approximation for
the spins correlations, $\langle s_{t}s_{0}\rangle=V_{t}$ . Then by iterating
this procedure, one can improve these approximations. Notice, however, that
once the noise correlation features a slow power-law decay, the spin
correlations exhibit a similar behavior because the spin dynamics are dictated
by the noise. Thus the spin correlations in a cycle of length $L$,

$c_{L}(t)=\left\langle s^{(0)}(0)s^{(0)}(t)\right\rangle_{L}$ (18)

(where averaging is only over cycles of length $L)$ decays as a power law
$t^{1-\alpha}$ within the range $L\ll t/\tau_{1}\ll L_{\max}$. In Fig. 5 we
demonstrate that this is indeed the behavior of $c_{L}(t)$ for $\alpha=3/2$
where correlations decay as $1/\sqrt{t}$. A full analysis of the behavior of
$c_{L}(t)$ in the limits $t\gg L$ and $|t|<L$ can be found in Appendices B &
C.

Figure 5: The long time asymptotic behavior of the spin correlation function
$c_{L}(t)$ for several cycles lengths, showing a power law decay. The
numerical results are obtained from 11000 realization of a system with
$2^{14}$ spins, $L_{\max}=2500$, $L_{\min}=1$, and $\alpha=\gamma=3/2$.

Consider now the case where the system is subjected to a constant magnetic
field, $h$. From Eq. (9) it follows that this amounts to a shift of the noise
average $\langle\zeta_{t}\rangle$ from zero to a finite value $h$, and
repeating the calculation that led to (10) now gives

$p(t)-\frac{1+m}{2}=V_{t}(h)\left[p(0)-\frac{1+m}{2}\right].$ (19)

Here $V_{t}(h)$ is given by an integral of the form (11) but with integration
range shifts such that, $-1<\zeta_{\nu}+h<1$, and $m$ is the magnetization
given by

$m=\frac{\mbox{erfc}\mu_{-}-\mbox{erfc}\mu_{+}}{\mbox{erfc}\mu_{-}+\mbox{erfc}\mu_{+}}$
(20)

where erfc$(x)$ is the complementary error function, and

$\mu_{\pm}=\frac{1\pm h}{\gamma\sqrt{2}\sqrt{1-\epsilon^{2}}}.$ (21)

In the strong coupling limit, $\gamma\gg 1$, and a weak magnetic field,
$h<\gamma$, one obtains that the magnetization saturates to

$M(t)\xrightarrow[t\to\infty]{}m\simeq\sqrt{\frac{2}{\pi}}\frac{h}{\gamma}.$
(22)

Notice that the shift dynamics within each cycle imply that the time scale for
polarizing the spins within a cycle of length $L$ is $\tau_{L}=\tau_{1}L$,
because the magnetic field affects only the connecting spins to other cycles.

A finite average magnetization changes noise properties. In particular, it
reduces the phase space for spin fluctuations. Therefore a transient noise
contribution, proportional to the square of the magnetization, is generated at
the turn ”ON” or turn ”OFF” of the external magnetic field. This contribution
can be neglected when the external magnetic field is small or when the
considered time interval is sufficiently far after the point where the
magnetic field has been turned off.

In principle, the above description applies only in the limit of large systems
$N\gg 1$, and strong coupling $\gamma\gg 1$, where relaxation is rapid. To
obtain slow relaxation, i.e., large $\tau_{1}$, we will consider situations
where $\gamma$ is of order one. However, $\gamma$ should be large enough to
avoid a situation where the cycles become, effectively, decoupled from each
other so that their magnetization freezes out. This limit is where the system
behavior is analogous to the glass phase of spin glasses. To avoid this phase,
$\gamma$ should be above 1.3. This value is obtained from numerical
simulations of the system. For cycles of equal length, $L=1$, a similar value
is obtained by requiring that the self-consistent equation (15) has a
solution.

##  V A Phenomenological Model

The mean field approach presented above suggests that the magnetization
$M_{L}(t)$ of each cycle (4) is an independent quantity which behaves
similarly to a noisy capacitor. Namely, the magnetization of a cycle of length
$L$ satisfies a Langevin equation (i.e. the Ornstein Uhlenbeck process):

$\frac{\partial M_{L}(t)}{\partial
t}=-\frac{M_{L}(t)-\tilde{h}(t)}{\tau_{L}}+\xi_{L}(t),$ (23)

where $\tilde{h}$ and $\xi_{L}(t)$ are the effective magnetic field and the
applied noise, respectively. In the strong coupling limit, $\gamma\gg 1$, and
weak magnetic field, $|h|\ll 1$, the relation between $\tilde{h}$ and the
magnetic field $h$ is given by equation (22):
$\tilde{h}(t)=\sqrt{2/\pi}h(t)/\gamma$. The noise, $\xi_{L}(t)$, is
approximated by a random Gaussian noise with zero mean and short range
correlations

$\langle\xi_{L}(t)\xi_{L}(t^{\prime})\rangle=\frac{2\sigma_{L}^{2}}{\tau_{L}}\delta(t-t^{\prime})$
(24)

where the constant $\sigma_{L}^{2}$ is determined such that the variance of
the magnetization obtained from the mean field approach equals that obtained
from the solution of the Langevin Eq. (23). The calculation of $\sigma_{L}$
can be found in Appendix D, see Eq. (86).

Taking the cycles to be independent, the distribution,
$P(M,t)=\langle\delta[M-M(t)]\rangle$, of the total magnetization (5), is
normal by the central limit theorem:

$P(M,t)=\frac{1}{\sqrt{2\pi}\sigma(t)}\exp\left[-\frac{[M-\overline{M}(t)]^{2}}{2\sigma^{2}(t)}\right].$
(25)

The mean magnetization at time $t$, $\overline{M}(t)$, and the variance,
$\sigma^{2}(t)$, can be calculated using Eq. (23).

Assuming the external magnetic field to be constant, $\tilde{h}$, within the
time range $(-t_{w},0)$ and zero otherwise, the solution of Eq. (23) for the
magnetization of a cycle of length $L$ is

$M_{L}(t)=\overline{M}_{L}(t)+\delta M_{L}(t),$ (26)

where for $t>0$,

$\overline{M}_{L}(t)=\tilde{h}\left[1-\exp\left(-\frac{t_{w}}{\tau_{L}}\right)\right]\exp\left(-\frac{t}{\tau_{L}}\right)$
(27)

is the average magnetization, while the fluctuating component is

$\delta
M_{L}(t)=\int_{-t_{*}}^{t}dt^{\prime}\xi_{L}(t^{\prime})\exp\left(-\frac{t-t^{\prime}}{\tau_{L}}\right).$
(28)

Here $-t_{*}$ is the initial time where the dynamics started.

The total average magnetization of the system is now obtained by averaging Eq.
(27) over the distribution of cycles length (1):

$\overline{M}(t)=\left\langle\overline{M}_{L}\right\rangle_{L}.$ (29)

This average can be expressed in terms of the exponential integral function,
$E_{n}(x),$

$\displaystyle\overline{M}(t)\\!$
$\displaystyle=\\!\nu\tilde{h}\left\\{\frac{1}{L_{\max}^{\alpha-1}}\\!\left[E_{2-\alpha}\\!\left(\frac{t}{\tau_{\max}}\right)\\!-\\!E_{2-\alpha}\left(\frac{t\\!+\\!t_{w}}{\tau_{\max}}\right)\right],\right.$
$\displaystyle-\left.\frac{1}{L_{\min}^{\alpha-1}}\\!\left[E_{2-\alpha}\\!\left(\frac{t}{\tau_{\min}}\right)\\!-\\!E_{2-\alpha}\left(\frac{t\\!+\\!t_{w}}{\tau_{\min}}\right)\right]\right\\}$
(30)

where $\nu$ is the normalization of the cycles length distribution function
(1), $\tau_{\max}=\tau_{1}L_{\max}$, and $\tau_{\min}=\tau_{1}L_{\min}$.

The curves that describe the average magnetization, $\overline{M}(t)$, in a
system with $2^{14}$ spins, $\alpha=\gamma=3/2$, $L_{\min}=1$,
$L_{\max}=2500$, and $h=0.8$, are depicted in Fig. 6 for various values of
waiting times $t_{w}$. The colored lines show the results of numerical
simulations, while the dashed lines are the analytic results of Eq. (30)
obtained when substituting $\tilde{h}=0.57$ and $\tau_{1}=1.9$. These values
are slightly different from those obtained from the large $\gamma$ asymptotic
expressions which are $\tau_{1}\simeq 1.68$ and $\tilde{h}\simeq 0.43$,
because $\gamma$ is not large and the external magnetic field, $h$, is not
weak.

Figure 6: A Log-log plot of the average magnetization of the RCCN model for
various waiting times. The colored lines are the results of simulation with
parameters as in Fig. 3, while the dashed lines are the curves obtained from
formula (30). (Only 5 values of waiting times, $t_{w}=20,40,160,640$, and
3000, are presented for clarity).

The variance of the fluctuations of the magnetization is obtained by squaring
Eq. (28) and averaging both over the noise $\zeta_{L}(t)$ and over the cycles
length distribution, while taking into account that we average over $N$
cycles. By neglecting the effect of magnetic field on the noise variance and
setting $t_{*}=0$, we obtain:

$\sigma^{2}(t)=\frac{1}{N}\left\langle\sigma_{L}^{2}\left[1-\exp\left(\frac{-2t}{\tau_{L}}\right)\right]\right\rangle_{L}.$
(31)

For our purpose in the next section we shall need only the saturated value of
this variance (at long time):

$\overline{\sigma}^{2}=\frac{\langle\sigma^{2}_{L}\rangle}{N}.$ (32)

(see Eq. (87) in the Appendix. D).

##  VI The survival probability

The typical behavior of a single realization of the magnetization after
turning the magnetic field off is illustrated in Fig. 7. It exhibits strong
fluctuations on top of slow decay. Due to these fluctuations, the
magnetization will cross the zero line base downwards and upwards in an
alternating manner. Our goal is to calculate the survival probability, $S(t)$
(or 1$-$cdf), i.e., the probability that the magnetization did not cross the
zero base line up to time $t$. Its relation to the distribution of times that
the magnetization first passes the zero line value, $\rho(t)$, (i.e. the
distribution of $T_{1}$ in Fig. 7) is described in Eq. (6) and can be also
written as

$\rho(t)=-\frac{\partial S(t)}{\partial t}.$ (33)

Figure 7: An illustration of the typical fluctuating behavior of the
magnetization (near the zero baseline) as it decays from a positive initial
value.

This problem is different from the crossing problem in the Ornstein-Uhlenbeck
process[44] because it involves strong temporal correlation of the
magnetization . Following [45] we define $P_{<}(t)$ to be the probability that
$M(t)<0$ at time $t$. From Eq. (25) it follows that

$P_{<}(t)=\int_{-\infty}^{0}dMP(M,t)=\frac{1}{2}\mbox{erfc}\left[\frac{\overline{M}(t)}{\sqrt{2}\sigma(t)}\right].$
(34)

To calculate the survival probability, we also need the distribution
functions, $\psi_{\pm}(\Delta t)$, of the time, $\Delta t,$ that the
magnetization returns to zero after it crossed from above and stayed negative
(“$-$”) or crossed from below and stayed positive (“$+$”). In terms of Fig. 7,
$\psi_{-}(\Delta t)$ is the distribution of the even time intervals $T_{2}$,
$T_{4}$, $T_{6}$ etc., while ,$\psi_{+}(\Delta t)$ is the distribution of the
odd time intervals $T_{3}$, $T_{5}$, $T_{7}$ etc. (Notice that the
distribution function for the first crossing at $T_{1}$, is special). The main
assumption that we need here is that all these time intervals are independent
random variables described by a distribution that changes adiabatically in
time. This assumption is based on our numerical simulations, which show that
the average magnetization changes very slowly compared to the typical time
intervals, $T_{k}$.

For a system with discrete-time evolution, the precise value of the crossing
time may be defined by interpolation, and the distribution functions
$\psi_{\pm}(\Delta t)$ should be averaged over the distribution of the
trajectory overshoot near the crossing point. Here we avoid this complexity by
setting our time resolution to be larger than a single time step (the typical
return time is of order 50-time steps); hence the magnetization of a single
realization is regarded as a continuous function.

Now, let us define $p_{2k-1}(t)$ to be the probability that a trajectory of
$M(t)$ starting at $M(0)>0$ ends below the zero baseline at time $t$ after
$2k-1$ crossing events (where $k\geq 1$). Then:

$P_{<}(t)=\sum_{k=1}^{\infty}p_{2k-1}(t)$ (35)

It would be instructive to examine $p_{2k-1}(t)$. If the first crossing event
took place at time $t_{1}<t$ and there were no other crossing events until
time $t$, then

$p_{1}(t)=\int_{0}^{t}dt_{1}\rho(t_{1})Q(t-t_{1}),$ (36)

where $Q(t)$ is the probability that $M(t^{\prime})$ remains negative for time
$t$. This is the probability of no crossing until time $t$, therefore it can
be written in terms of the return probability, i.e.

$Q(t)=1-\int_{0}^{t}dt^{\prime}\psi_{-}(t^{\prime}).$ (37)

Consider now the next function $p_{3}(t)$. Assuming the downward crossing took
place at $t_{1}$, the upward crossing at $t_{2}$ and the next downward
crossing at $t_{3}$ so that $0\leq t_{1}\leq t_{2}\leq t_{3}\leq t$ we have,

$\displaystyle
p_{3}(t)=\int_{0}^{t}dt_{1}\rho(t_{1})\int_{t_{1}}^{t}dt_{2}\psi_{-}(t_{2}-t_{1})$
(38)
$\displaystyle\times\int_{t_{2}}^{t}dt_{3}\psi_{+}(t_{3}-t_{2})Q(t-t_{3}).$

Similar convolution integrals describe $p_{2k-1}(t)$ with $k>2$. Taking the
Laplace transform of these formulas we have

$\hat{p}_{2k-1}(s)=\hat{\rho}(s)[\hat{\psi}_{-}(s)\hat{\psi}_{+}(s)]^{k-1}\hat{Q}(s),$
(39)

were $\hat{f}(s)$ denotes the Laplace transform of $f(t)$.

From the above definitions it follows that

$P_{<}(t)=\sum_{k=1}^{\infty}p_{2k-1}(t),$ (40)

and taking the Laplace transform of equation (35), using (39), and solving for
$\hat{\rho}(s),$ we obtain:

$\hat{\rho}(s)=s\hat{P}_{<}(s)\frac{1-\hat{\psi}_{-}(s)\hat{\psi}_{+}(s)}{1-\hat{\psi}_{-}(s)}.$
(41)

Finally, from the above equation and the Laplace transform of Eq. (33), with
$S(0)=1,$ we obtain

$\hat{S}(s)=\frac{1}{s}-\hat{P}_{<}(s)\frac{1-\hat{\psi}_{-}(s)\hat{\psi}_{+}(s)}{1-\hat{\psi}_{-}(s)}.$
(42)

In order to proceed, one needs the functional form of $\psi_{\pm}(t)$. These
functions are not known and difficult to compute. Here we assume that each of
these functions is characterized by a single time-scale and has an exponential
form:

$\psi_{\pm}(t)=\frac{1}{\tau_{\pm}}\exp\left(-\frac{t}{\tau_{\pm}}\right),$
(43)

where $\tau_{\pm}$ are the decay time scales. Taking the Laplace transform of
these function and substituting them in Eq. (42) we obtain:

$\hat{S}(s)=\frac{1}{s}-\hat{P}_{<}(s)-\frac{\tau_{+}}{\tau_{-}}\frac{1}{1+\tau_{+}s}\hat{P}_{<}(s).$
(44)

In order to identify the ratio $\tau_{+}/\tau_{-},$ we multiply the above
equation by $s$ and take the limit $s\to 0$. Then using the following property
of Laplace transform: $\lim_{s\to 0}s\hat{f}(s)=\lim_{t\to\infty}f(t)$, and
the assumption of adiabatic evolution, namely that the magnetization
$\overline{M}(t)$ and the standard deviation $\sigma(t)$ can be assumed to be
constants over large time intervals (compared to $\tau_{\pm}$), we obtain:

$\frac{\tau_{+}}{\tau_{-}}=\frac{1}{P_{<}(t)}-1$ (45)

Thus, the inverse Laplace transform of Eq. (44), leads to

$S(t)=1-P_{<}(t)-\left(\frac{1}{P_{<}(t)}-1\right)\int_{0}^{t}\frac{dt^{\prime}}{\tau_{+}}e^{-t^{\prime}/\tau_{+}}P_{<}(t-t^{\prime}).$
(46)

The decay time $\tau_{+}$ is unknown. From the adiabatic assumption it is
clear that it should reduce slowly in time because as the average
magnetization becomes smaller the probability for $M(t)$ to stay positive
reduces. Moreover, ageing suggests that $\tau_{+}$ also has a weak dependence
on the waiting time, $t_{w}$ which sets the value of the average
magnetization, see Appendix E. The dashed lines in Fig. 3 (corresponding to
$t_{w}=20,40,160,640$ and 3000) are obtained from formula (46) with
logarithmic dependence on the waiting time: $\tau_{+}=11\log(t_{w}/2)$ (see
Appendix E). The other parameters are $\tau_{1}=1.9$, $\tilde{h}=0.57$ (the
same as for the average magnetization shown in Fig. 6), and
$\sigma(t)\simeq\overline{\sigma}=0.047$ (see Appendix D).

##  VII Conclusion

To conclude, we constructed a minimal toy model that reproduces the
statistical behavior of the division time of a E. coli cell that undergoes a
period of acute stress in which its growth is arrested. This construction
rests on the following central assumptions: First, the cellular network may be
considered random when taking the cell far from its adaptive regime. Second,
one can describe this network using Boolean variables that account for each
cell component’s active or inactive states. Third, the statistical properties
of the dynamics of the system are captured by discrete-time evolution. Fourth,
the interactions among components can be approximated by two-body
interactions. Namely, the state of a spin depends (nonlinearly) only on sums
of the form $\sum_{j}J_{ij}s_{j}$ and not, e.g., terms of the form
$\sum_{j,k}J_{ijk}s_{j}s_{k}$ that represent three-body interactions. With
these assumptions, the simplest model we could construct that reproduces the
experimental results is the RCCN model illustrated in Fig. 2.

Clearly, none of these assumptions hold on a microscopic level description:
The cellular network is not random, as it evolved during billions of years of
evolution; the concentration of cell components and the time evolution are
continuous; and reactions within the cell usually involve few-body
interactions. Nevertheless, the agreement between the experimental data and
the numerical and analytical solution of the model suggests that the RCCN
model provides an effective statistical description of the cell. In other
words, one expects that a realistic microscopic model of the cellular network
can be reduced to the RCCN model in the proper regime of parameters. As such,
this model forms a starting point on which a more refined description can be
built by taking into account particular features of the cellular network.

Moreover, our study shows that the random statistical behavior of the cells
holds valuable information about the cellular network structure. For instance,
comparing the solution of the RCCN model with the experimental data allows one
to extract information about the cycles length distribution, e.g., the power
$\alpha$ in Eq. (1), as well as the scale of the longest cellular cycle,
$L_{\max}$ of the cell.

In order to see how $\alpha$ is reflected in the survival probability,
consider the intermediate asymptotic behavior of $S(t)$ within the range
$\tau_{\min}\ll t\ll\tau_{\max}$ which can be evaluated assuming $\tau_{+}\ll
t$. Within this range the upper limit in the integral (46) can be extended to
infinity, while the function $P_{<}(t-t^{\prime})$ can be expanded in powers
of $t^{\prime}$. The zeroth order term of this expansion cancels the first
term on the right hand side of Eq. (46), while the first order term yields,

$S(t)\approx\left\\{\begin{array}[]{lr}b_{+}\left(\frac{\tau_{\max}}{t}\right)^{\alpha+1}&\tau_{\min}\ll
t_{w}\ll t\ll\tau_{\max}\\\
b_{-}\left(\frac{\tau_{\max}}{t}\right)^{\alpha}&\tau_{\min}\ll
t\ll\tau_{\max}<t_{w}\end{array}\right.,$ (47)

with

$b_{\pm}=(\alpha-1)\Gamma(\alpha)\frac{\tilde{h}}{\sqrt{2\pi}\overline{\sigma}}\left(\frac{L_{\min}}{L_{\max}}\right)^{\alpha-1}\frac{\tau_{+}}{\tau_{\max}}\left(\frac{\alpha
t_{w}}{\tau_{\max}}\right)^{\mu_{\pm}}$ (48)

where $\mu_{+}=1$, and $\mu_{-}=0$.

Eq. (47) shows that the distribution of the cycles length (1) determines the
intermediate asymptotic behavior of the survival probability: The power-law
decay changes from $t^{-\alpha-1}$ for short waiting time to $t^{-\alpha}$ at
saturation, i.e. when $t_{w}\geq\tau_{\max}$. Notice that this result is
independent of any other parameter of the system. Thus, the behavior of the
survival probability at saturation reflects the distribution of the lengths of
the cycles [11].

From a theoretical viewpoint, the RCCN model is a complicated model, and the
phenomenological description that we have presented here provides only a basic
description. For example, we did not take the effect of magnetization on the
noise level, the precise dependence of $\tau_{+}$ on time, and we also did not
account for possible cross-correlations between different cycles whose
cumulative effect, despite being very small, may become significant.

The theoretical framework provided by the RCCN model opens the possibility to
study many other universal features of the system such as correlations between
cycles; the effect of stress focused on a single node or a limited number of
them; and the sensitivity of the dynamics to a small change in the initial
conditions. These questions are left for future study.

##  VIII Appendix A: Ergodicity

Eq. (8) is straightforwardly obtained by ensemble averaging. In the first part
of this Appendix, we show that the same equation can also be derived by
averaging over initial conditions of the spin configurations. However, this
derivation relies on the assumption that different cycles are uncorrelated. In
the second part of the Appendix, we show that the correlation between cycles
decay as $1/\sqrt{N}$ where $N$ is the number of cycles.

In what follows we use overbar to denote averaging over initial conditions and
over the cycles. Using definition (7), we have

$\overline{\zeta_{i}(t)\zeta_{i}(t^{\prime})}=\frac{1}{N}\sum_{i}\sum_{k\neq
i}\sum_{k^{\prime}\neq
i}J_{ik}J_{ik^{\prime}}\overline{s_{k}(t)s_{k^{\prime}}(t^{\prime})}.$ (49)

Hereinafter we suppress the upper index that denote the connecting spins in
each cycles. Assuming the cycles to be uncorrelated implies that only terms
with $k=k^{\prime}$ contribute to the above sum. This diagonal approximation
yields

$\overline{\zeta_{i}(t)\zeta_{i}(t^{\prime})}=\frac{1}{N}\sum_{i}\sum_{k\neq
i}J_{ik}^{2}\overline{s_{k}(t)s_{k}(t^{\prime})}$ (50)

and by the law of large numbers, in the limit $N\to\infty$, we obtain

$\overline{\zeta_{i}(t)\zeta_{i}(t^{\prime})}=\gamma^{2}\overline{s_{k}(t)s_{k}(t^{\prime})}.$
(51)

Comparing this formula with Eq.(8) shows that averaging over initial
conditions is equivalent to ensemble averaging provided
$\overline{s_{k}^{(0)}(t)s_{k^{\prime}}^{(0)}(t^{\prime})}\to 0$ when
$N\to\infty$ and $k\neq k^{\prime}$. To validate this property consider the
following quantity:

$\Phi_{ij}=\left\langle
J_{ij}\overline{s_{i}(t)s_{j}(t)}\right\rangle,~{}i\neq j.$ (52)

Expressing $s_{i}(t)$ in terms of its previous time step, using Eqs. (2) we
have

$\Phi_{ij}=\left\langle
J_{ij}\overline{\mbox{sign}[s_{i}(t-L_{i})+\zeta_{i}(t-1)]s_{j}(t)}\right\rangle,$
(53)

where $\zeta_{i}(t-1)$ is the input signal to the $i$-th cycle . Let us
separate from this signal the contribution that comes from the $j$-th spin,
$\zeta_{i}(t)=J_{ij}s_{j}(t)+\tilde{\zeta}_{i}(t)$, and expand sign function
to linear order in $J_{ij}$ (this expansion is justified because $J_{ij}$ is
of order $1/\sqrt{N}$ and we consider the limit $N\to\infty$). Thus

$\displaystyle\Phi_{ij}$ $\displaystyle=\left\langle
J_{ij}\overline{\mbox{sign}[s_{i}(t-L_{i})+J_{ij}s_{j}(t-1)+\tilde{\zeta}_{i}(t-1)]s_{j}(t)}\right\rangle$
$\displaystyle=\left\langle
J_{ij}\overline{\mbox{sign}[s_{i}(t-L_{i})+\tilde{\zeta}_{i}(t-1)]s_{j}(t)}\right\rangle$
$\displaystyle+\left\langle
J_{ij}^{2}\overline{s_{j}(t-1)s_{j}(t)2\delta[s_{i}(t-L_{i})+\tilde{\zeta}_{i}(t-1)]}\right\rangle$

Now let us assume, self consistently, that the spins are uncorrelated. Then
$\tilde{\zeta}_{i}(t)$ is a Gaussian noise, independent of $s_{i}$ and
$s_{j}$, with zero mean and with variance $\gamma^{2}$; hence, taking first
the mean over the noise only,
$\overline{\mbox{sign}[s_{i}(t)+\tilde{\zeta}_{i}]}=s_{i}(t)\overline{\mbox{sign}[1+\tilde{\zeta}_{i}]}$
and the term in the middle line of Eq. (VIII) vanishes upon averaging.
Similarly $\overline{\delta[s_{i}(t)+\tilde{\zeta}_{i}]}=P[-s_{i}(t)]$ where
$P(\tilde{\zeta}_{i})$ is the noise distribution function. Taking into account
that the maximal value of this distribution function is
$1/(\sqrt{2\pi}\gamma)$ and that the maximal value of $s_{j}(t)$ is $1$, we
obtain

$\Phi_{ij}\leq\left\langle
J_{ij}^{2}\max[P(\tilde{\zeta})]\right\rangle=\frac{\gamma}{\sqrt{2\pi}N}$
(55)

Noting that the spatial correlations is fully determined by $J_{ij}$, This
result together with the definition of (52) implies that

$\overline{s_{i}(t)s_{j}(t)}\sim\frac{J_{ij}+J_{ji}}{\gamma}\sim\frac{1}{\sqrt{N}}.$
(56)

##  IX Appendix B: Derivation of Eq. (10)

Consider the mean field approximation for the RCCN model where all cycles are
of length $L=1$. The spin dynamical equations are:

$s_{t+1}=\mbox{sign}(s_{t}+\zeta_{t})$ (57)

where $\zeta_{t}$ is the input noise. Let us look first on the case where the
noise is uncorrelated. Then the probability, $p(t)$ that $s_{t}=1$ satisfies
the equation:

$p(t+1)=p(t)P(\zeta_{t}>-1)+[1-p(t)]P(\zeta_{t}>1)$ (58)

Here the first term describes the probability that $s_{t}=1$ multiplied by the
probability that the noise cannot flip it, i.e., the probability that
$\zeta_{t}>-1$. The second term describes the case where $s_{t}=-1$ but the
noise is sufficiently large to flip it. Next, we use the following general
properties of the noise probabilities:

$P(\zeta_{t}>-1)=q+P(\zeta_{t}>1)$ (59a) where $q=P(-1<\zeta_{t}<1),$ (59b)
and $P(\zeta_{t}>1)+P(\zeta_{t}<-1)+q=2P(\zeta_{t}>1)+q=1,$ (59c)

where here we use the property that the noise distribution is an even
function, i.e. $P(\zeta_{t}>1)=P(\zeta_{t}<-1)$. Substituting these equations
in (58) gives

$p(t+1)-\frac{1}{2}=q\left[p(t)-\frac{1}{2}\right]$ (60)

From here Eq. (10) directly follows with $V_{t}=q^{t}$.

Consider now the noise be correlated and let us calculate the probability that
after two steps $s_{2}=1$. This probability is given by the following sum:

$p(2)=P(\zeta_{1}>1)+P(-1<\zeta_{1}<1,\zeta_{0}>1)+P(-1<\zeta_{1}<1,-1<\zeta_{0}<1)p(0).$
(61)

The first term on the right-hand side of this equation is associated with the
situation that the noise at the last step is sufficiently large to ensure that
the spin will be one independent of its state. The second term describes the
case where the noise in the last step is within the range where no spin flip
occur (irrespective of the initial spin state), but the noise in the previous
step is large enough to ensure that the spin is in ”one” state. Finally, the
last term is the probability that the spin was initially at one state, and the
noise in both time steps doesn’t flip the spin. These probabilities refer to
disjoint sets of events and therefore add up.

Now, using relations similar to (59) allows one to rewrite the above equation
in the form:

$p(2)=\frac{1}{2}\left[1-P(-1<\zeta_{1}<1)\right]+\frac{1}{2}[P(-1<\zeta_{1}<1)-V_{2}]+V_{2}p(0),$
(62)

where

$V_{2}=P(-1<\zeta_{1}<1,-1<\zeta_{0}<1).$ (63)

Thus

$p(2)-\frac{1}{2}=V_{2}\left[p(0)-\frac{1}{2}\right].$ (64)

This derivation can be generalized to any number of time steps giving Eq.
(10).

##  X Appendix C: Noise in the RCCN model

In this Appendix, we study the fluctuation of the cycles magnetization. To
this end, we shall characterize the correlation function (18) of the
connecting spin in a cycle of length $L$, and use it to calculate the
magnetization variance in order to identify $\sigma_{L}^{2}$, and its average
over the distribution of the lengths of the cycles, $\overline{\sigma}^{2}$.

The essential quantity needed for our analysis is the conditional probability
of a connecting spin and the noise acting on the same spin at some different
time, $P(\zeta_{t}|s_{0})$, where the subscript denotes the time step. It is
easy to calculate this conditional probability for $t=-1$ because the time
evolution $s_{0}=\mbox{sign}(s_{-1}+\zeta_{-1})$ allows one to calculate the
conditional probability $P(s_{0}|\zeta_{-1})$. Then using Bayes’ theorem, we
have

$P(\zeta_{t-1}|s_{0})=P(s_{0}|\zeta_{-1})\frac{P_{0}(\zeta_{-1})}{P(s_{0})}$
(65)

where $P(s_{0})=1/2$ is the unconditional probability of the spin (in the
absence of magnetic field), and $P_{0}(\zeta)$ is the probability distribution
of the noise at some arbitrary time, which is normal distribution with zero
mean and variance $\gamma^{2}$. From here we obtain that for $s_{0}=1$:

$P(\zeta_{-1}|1)=\left\\{\begin{array}[]{cc}2P_{0}(\zeta_{-1})&\zeta_{-1}>1\\\
P_{0}(\zeta_{-1})&|\zeta_{-1}|<1\\\ 0&\zeta_{-1}<-1\end{array}\right.$ (66)

The conditional probability is symmetric with respect to a change sign of both
$s_{0}$ and $\zeta_{-1}$ hence $P(\zeta_{-1}|-1)=P(-\zeta_{-1}|1)$.

Now we can calculate the conditional probability at any other time by the
integral

$P(\zeta_{t}|s_{0})=\int
d\zeta_{-1}P(\zeta_{t}|\zeta_{-1},s_{0})P(\zeta_{-1}|s_{0}).$ (67)

Here $P(\zeta_{t}|\zeta_{-1},s_{0})$ is the conditional probability of the
noise at time $t$ given the noise value at time $t=-1,$ and the state of the
spin, $s_{0}$. However, the effect of a single spin on the noise is negligible
because the noise results from a very large number of spins; hence

$\displaystyle P(\zeta_{t}|\zeta_{-1},s_{0})$ $\displaystyle=$ $\displaystyle
P(\zeta_{t}|\zeta_{-1})$ (68) $\displaystyle=$
$\displaystyle\frac{\exp\left(-\frac{(\zeta_{t}-\mu_{t+1}\zeta_{-1})^{2}}{2\gamma^{2}(1-\mu_{t+1}^{2})}\right)}{\gamma\sqrt{2\pi(1-\mu_{t+1}^{2})}}$

where $\mu_{t+1}=\langle\zeta_{t}\zeta_{-1}\rangle/\gamma^{2}$ is the noise
correlation function normalized by its variance. Substituting (66) and (68) in
(67) yields:

$P(\zeta_{t}|s_{0})=\left[1+s_{0}\chi(\zeta_{t})\right]P_{0}(\zeta_{t})$ (69)

where

$\chi(\zeta)=\frac{1}{2}\sum_{\pm}\pm\mbox{erf}\left(\frac{1\pm\zeta\mu_{t+1}}{\gamma\sqrt{2-2\mu_{t+1}^{2}}}\right).$
(70)

We turn now to calculate the autocorrelation function of the connecting spin
in a cycle of length $L$ at time $t=L$, $c_{L}(L)$, in the limit of zero
magnetic field. Setting $t=L-1$ in Eq. (9), multiplying it by $s_{0}$, and
taking the average we obtain

$c_{L}(L)=\langle s_{0}s_{L}\rangle=\langle
s_{0}\mbox{sign}(s_{0}+\zeta_{L-1})\rangle.$ (71)

This average can be evaluated using (69) giving

$c_{L}(L)=\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)+2\int_{1}^{\infty}d\zeta\chi(\zeta)P_{0}(\zeta).$
(72)

This is a general function of $\mu_{L}$ and $\gamma$, but for long cycles one
can assume that $\mu_{L}\ll 1$, and the above integral may be expanded up to
linear order in this parameter, giving

$c_{L}(L)\simeq\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)+q\mu_{L},$ (73)

where

$q=\frac{2}{\pi}\exp\left(-\frac{1}{\gamma^{2}}\right).$ (74)

Consider, now the same correlation function but at time $t=2L$, i.e.,
$c_{L}(2L).$ Setting $t=2L-1$ in Eq. (9), multiplying by $s_{0}$ , and
expressing $s_{L}$ in terms of $s_{0}$ using Eq.(9) once again we obtain:

$c_{L}(2L)=\left\langle
s_{0}~{}\mbox{sign}[~{}\mbox{sign}(s_{0}+\zeta_{L-1})+\zeta_{2L-1}]\right\rangle.$
(75)

To evaluate this quantity, we need the joint distribution of the noise at two
points in time conditioned by the state of the spin $s_{0}$:

$P(\zeta_{t},\zeta_{t^{\prime}}|s_{0})=\int d\zeta_{-1}\
P_{0}(\zeta_{t},\zeta_{t^{\prime}}|\zeta_{-1})P(\zeta_{-1}|s_{0})$ (76)

To leading order in the noise correlations, $\mu_{t}$, This integral gives:

$\frac{P_{0}(\zeta_{t},\zeta_{t^{\prime}}|1)}{P_{0}(\zeta_{t})P_{0}(\zeta_{t^{\prime}})}\simeq
1+\zeta_{t}\zeta_{t^{\prime}}\mu_{t-t^{\prime}}+\frac{\sqrt{q}}{\gamma}(\zeta_{t}\mu_{t}+\zeta_{t^{\prime}}\mu_{t^{\prime}}).$
(77)

The average on the right hand side of (75) can be evaluated by expressing the
outer sign function in terms of its Fourier integral,

$\mbox{sign}(x)=2\mbox{Re}\int_{0}^{\infty}\frac{d\eta}{i\pi\eta}\exp(i\eta
x).$ (78)

Preforming the average yields

$c_{L}(2L)\simeq\mbox{erf}^{2}\left(\\!\frac{1}{\sqrt{2}\gamma}\\!\right)+q\left[\mu_{L}\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)+\mu_{2L}\right].$
(79)

Notice that the contribution from the noise autocorrelation, i.e the quadratic
term on the right-hand side of Eq. (77), vanishes. Thus one can interpret the
above formula in the following way: The first term is simply the probability
that the spin does not flip by the noise after two events where its state may
have been changed (at time $L$ and time $2L)$, assuming these events to be
independent. The two other contributions represent corrections to this result
due to correlations of the noise with the spin. The first is the product of
the probabilities that noise does not flip the spin at one time event
multiplied by the excess probability for the spin to remain in its state due
to correlations with the spin after time $L$. The second contribution comes
from correlations of the noise and the spin after time $2L$.

A similar calculation for $n$ time steps (of length $L$), in the asymptotic
limit $L\gg 1$ gives

$c_{L}(nL)\simeq
q\sum_{j=1}^{n}\mbox{erf}^{n-j}\left(\\!\frac{1}{\sqrt{2}\gamma}\right)\mu_{jL}+\mbox{erf}^{n}\left(\\!\frac{1}{\sqrt{2}\gamma}\\!\right).$
(80)

Assuming this behavior is characteristic for any discrete time $t$ and not
only to the values $t=nL$ one can now solve this equation self consistently by
averaging $c_{L}(t)$ over the power law distribution of the cycles length (1).
This average should be proportional to the noise autocoronation $\mu_{t}$.

Consider, first, the average of the last term in the above formula (which
decays exponentially). In the asymptotic limit $\tau_{\min}\ll
t\ll\tau_{\max}$ it yields

$\left\langle\mbox{erf}^{\frac{t}{L}}\left(\frac{1}{\sqrt{2}\gamma}\right)\right\rangle_{L}=\nu\Gamma(\alpha-1)\left(\frac{\tau_{1}}{t}\right)^{\alpha-1}$
(81)

where in the approximation where correlations are taken only to first order
to, $\tau_{1}^{-1}\simeq-\log\mbox{[erf}(1/\sqrt{2}\gamma)]$. This term yield
an algebraic decay, and substituting the same form of decay, $\mu_{jL}\simeq
b(jL)^{1-\alpha}$, in Eq. (80), the sum can be expressed in terms of Lerch
transcendent which yields the following asymptotic formula:

$c_{L}(t)\simeq\tilde{b}\left(\frac{1}{t+L}\right)^{\alpha-1}\\!\\!+\exp\left(-\frac{t}{\tau_{1}L}\right)$
(82)

where

$\tilde{b}=\frac{bq{}}{1-\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)}.$
(83)

This result shows that in the large time asymptotic limit, all cycles
correlation functions decay as a power law with the same power and amplitude
as illustrated in Fig. 5 and Fig. 8.

##  XI Appendix D: The magnetization noise

In this appendix, we calculate the magnetization noise variance,
$\overline{\sigma}^{2}$, at zero magnetic field. It is obtained by averaging
the variance of the magnetization noise, $\sigma_{L}^{2}$, of cycles of length
$L$ over the cycles length distribution (1).

Consider the magnetization of a cycle of length $L$. Taking into account that
the dynamics within the cycle is shift dynamics, we have

$M_{L}(t)=\frac{1}{L}\sum_{k=0}^{L-1}s^{(k)}(t)=\frac{1}{L}\sum_{t^{\prime}=0}^{L-1}s^{(0)}(t-t^{\prime}).$
(84)

Thus when the system reaches a stationary state, the variance of the
magnetization of a cycle of length $L$ is given by:

$\displaystyle\langle M_{L}^{2}\rangle$
$\displaystyle=\frac{1}{L^{2}}\sum_{t^{\prime\prime}=0}^{L-1}\sum_{t^{\prime}=0}^{L-1}\langle
s^{(0)}(t-t^{\prime\prime})s^{(0)}(t-t^{\prime})\rangle$
$\displaystyle=\frac{1}{L^{2}}\sum_{t^{\prime\prime}=0}^{L-1}\sum_{t^{\prime}=0}^{L-1}c_{L}(t^{\prime\prime}-t^{\prime})$
$\displaystyle=\frac{1}{L}+\frac{2}{L^{2}}\sum_{t=1}^{L-1}(L-t)c_{L}(t),$ (85)

where we took into account that $c_{L}(0)=1$.

The above sum depends on the behavior of $c_{L}(t)$ within the short time
regime $0<t<L$. The results obtained in Appendix C (see Eqs. (82) and (73))
indicate that within this range, the correlation length is essentially
constant. This conclusion is also supported by the numerical results shown in
Fig. 8, where we depict $c_{L}(t)$ at three values of cycle length for
$\alpha=3/2$. This figure shows modulations of the correlations on the scale
of $L$ with an approximate constant behavior within the range $1\leq t\leq
L-1$.

If we set $c_{L}(t)$ to be constant $k$, then

$\sigma_{L}^{2}\simeq\langle M_{L}^{2}\rangle\simeq\frac{1}{L}+k\frac{L-1}{L}$
(86)

To calculate the $\overline{\sigma}^{2}$ in Eq. (32), one should take the
average of $\sigma_{L}^{2}$ over the cycles length distribution. For
simplicity we assume $L_{\min}=1$ and extend the sum to infinity (as it
converges for $\alpha>1$). With these approximations we obtain

$\overline{\sigma}^{2}=\frac{1}{N}\langle\sigma_{L}^{2}\rangle=\frac{1}{N}\
\left[k+\frac{(1-k)\zeta(1+\alpha)}{\zeta(\alpha)}\right],$ (87)

where $N$ is the average number of cycles, while $\zeta(x)$ is the Riemann
zeta function.

Finally, it remains to obtain an estimate for the constant $k$. For this
purpose, let us calculate the correlation function at one time step,

$c_{L}(1)=\langle s_{0}s_{1}\rangle=\langle
s_{0}\mbox{sign}(s_{1-L}+\zeta_{1})\rangle.$ (88)

In order to calculate this average one needs the conditional probability
$P(\zeta_{1}|s_{0},s_{1-L})$, namely the distribution of the noise at time
$t=1$ for a given state of the spin at previous times $t=0$ and $t=1-L$. For
sufficiently large cycles, one expects that the correlations of the noise with
the state of the spin at the latter time can be neglected and to use the
approximation $P(\zeta_{1}|s_{0},s_{1-L})\simeq P(\zeta_{1}|s_{0})$. Then
using (69) to perform the average, and linearizing with respect to $\mu_{1}$
(although this is not strictly justified because $\mu_{1}$ is of order one,
but here we are only interested in the approximate value) we obtain:

$c_{L}(1)=\left\langle
s_{0}\mbox{erf}\left(\frac{s_{1-L}}{\sqrt{2}\gamma}\right)\right\rangle+q\mu_{1}.$
(89)

Next, to perform the average over $s_{1-L}$ we use the conditional
probability:

$P(s_{1-L}|s_{0})=\left\\{\begin{array}[]{ll}\frac{1+u}{2}&\mbox{
for}~{}~{}s_{1-L}=s_{0}\\\ \frac{1-u}{2}&\mbox{
for}~{}~{}s_{1-L}=-s_{0}\end{array}\right.$ (90)

where $u=\langle s_{0}s_{1-L}\rangle=c_{L}(L-1)$. From the last two equations
we obtain:

$\displaystyle c_{L}(1)$ $\displaystyle\simeq
c_{L}(L-1)\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)+q\mu_{1}$ (91)

This formula connects the correlation function at one time step with that at
the time step just before the completion of a full cycle of $L$ steps. Setting
$k=c_{L}(1)\simeq c_{L}(L-1)$ and solving for $k$ we obtain.

$k\simeq\frac{q\mu_{1}}{1-\mbox{erf}\left(\frac{1}{\sqrt{2}\gamma}\right)}.$
(92)

Figure 8: A log-log plot of a connecting spin correlation function,
$c_{L}(t)$, for various values of cycles length. The dots are connected by
lines for clarity.

In our numerical study $\alpha=\gamma=3/2,$ which implies that $k\simeq
0.81\mu_{1}$. Taking also that average number of cycles is $N=408$, and
approximating $\mu_{1}\simeq 1$ we obtain from (87) that
$\overline{\sigma}=0.047$ which is the same as the numerical result.

##  XII Appendix E: The waiting time Dependence of $\tau_{+}$

This appendix explains the logarithmic dependence of the return time,
$\tau_{+},$ on the waiting time, $t_{w}$. For this purpose, we take note that
formula (46) is obtained by taking into account correlations in the relaxation
of the magnetizaion and the effect of the second term is manifested only at a
sufficiently long time, when the survival probability becomes small. Thus in
what follows, we estimate $\tau_{+}$ by focusing our attention on the regime:

$t\sim\tau_{\max},\quad\mbox{and}\quad\tau_{\min}\ll t_{w}<\tau_{\max}.$ (93)

The magnetization is a sum of two contributions: the mean value of the
magnetization and a fluctuating part. From Eq. (30) it follows that average
magnetization in the regime of interest (93) is given by

$\overline{M}\approx ct_{w}.$ (94)

where

$c\simeq(\alpha-1)\Gamma(\alpha,1)\frac{\tilde{h}}{\tau_{\max}}\left(\frac{L_{\min}}{L_{\max}}\right)^{\alpha-1},$
(95)

and $\Gamma(a,z)$ is the incomplete gamma function. This positive contribution
comes mainly from the long cycles which decay slowly and keep memory of the
time that the system was subjected to the magnetic field. It is therefore
clear that it increases with the waiting time, $t_{w}$, (as long as it is
shorter than the longest cycle).

On the other hand, the main contribution to the fluctuating part of the
magnetization comes from short cycles that quickly lose memory of the magnetic
field, and because there are many of them. The relaxation of these very short
orbits is approximately exponential as follows from (27).

Now consider a trajectory of the magnetization that crosses zero value
upwards. Since the average magnetization is positive, the fluctuation must be
negative with amplitude given by (94). Thus one expects it to have the form
$\delta M(t)\approx-\bar{M}\exp(-t/\tau_{\min})$, where time is measured from
the crossing point. Now, the typical time that the magnetization remains
positive is several times (say $n$) of the time that takes this fluctuation to
relax. Namely it is obtained from the condition $\delta M(\tau_{+}/q)=-\delta$
where $\delta$ is some small (positive) value of the magnetization. Solution
of this equation yields a logarithmic dependence on the return time,
$\tau_{+}=n\tau_{\min}\log(ct_{w}/\delta)$.

## References

  * Parisi [1993] G. Parisi, Statistical physics and biology, Physics World 6, 42 (1993).
  * Ozbudak  _et al._ [2004] E. M. Ozbudak, M. Thattai, H. N. Lim, B. I. Shraiman, and A. Van Oudenaarden, Multistability in the lactose utilization network of escherichia coli, Nature 427, 737 (2004).
  * Milo  _et al._ [2002] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, Network motifs: simple building blocks of complex networks, Science 298, 824 (2002).
  * Kauffman  _et al._ [2003] S. Kauffman, C. Peterson, B. Samuelsson, and C. Troein, Random boolean network models and the yeast transcriptional network, Proceedings of the National Academy of Sciences 100, 14796 (2003).
  * Li  _et al._ [2013] Z. Li, S. Bianco, Z. Zhang, and C. Tang, Generic properties of random gene regulatory networks, Quantitative biology 1, 253 (2013).
  * Edwards and Glass [2000] R. Edwards and L. Glass, Combinatorial explosion in model gene networks, Chaos: An Interdisciplinary Journal of Nonlinear Science 10, 691 (2000).
  * Stern  _et al._ [2014] M. Stern, H. Sompolinsky, and L. Abbott, Dynamics of random neural networks with bistable units, Physical Review E 90, 062710 (2014).
  * Stern  _et al._ [2007] S. Stern, T. Dror, E. Stolovicki, N. Brenner, and E. Braun, Genome-wide transcriptional plasticity underlies cellular adaptation to novel challenge, Molecular Systems Biology 3, 106 (2007).
  * Braun [2015] E. Braun, The unforeseen challenge: from genotype-to-phenotype in cell populations, Reports on Progress in Physics 78, 036602 (2015).
  * Schreier  _et al._ [2017] H. I. Schreier, Y. Soen, and N. Brenner, Exploratory adaptation in large random networks, Nature communications 8, 1 (2017).
  * Kaplan  _et al._ [2021] Y. Kaplan, S. Reich, E. Oster, S. Maoz, I. Levin-Reisman, I. Ronin, O. Gefen, O. Agam, and N. Q. Balaban, Observation of universal ageing dynamics in antibiotic persistence, Nature 600, 290 (2021).
  * Tosa and Pizer [1971] T. Tosa and L. I. Pizer, Effect of serine hydroxamate on the growth of escherichia coli, Journal of bacteriology 106, 966 (1971).
  * Ritort and Sollich [2003] F. Ritort and P. Sollich, Glassy dynamics of kinetically constrained models, Advances in physics 52, 219 (2003).
  * Struik [1977] L. C. E. Struik, Physical aging in amorphous polymers and other materials, (1977).
  * Hodge [1995] I. M. Hodge, Physical aging in polymer glasses, Science 267, 1945 (1995).
  * Hwa  _et al._ [2003] T. Hwa, E. Marinari, K. Sneppen, and L.-h. Tang, Localization of denaturation bubbles in random dna sequences, Proceedings of the National Academy of Sciences 100, 4411 (2003).
  * Matan  _et al._ [2002] K. Matan, R. B. Williams, T. A. Witten, and S. R. Nagel, Crumpling a thin sheet, Physical Review Letters 88, 076101 (2002).
  * Lahini  _et al._ [2017] Y. Lahini, O. Gottesman, A. Amir, and S. M. Rubinstein, Nonmonotonic aging and memory retention in disordered mechanical systems, Physical review letters 118, 085501 (2017).
  * Ghofraniha  _et al._ [2007a] N. Ghofraniha, C. Conti, and G. Ruocco, Aging of the nonlinear optical susceptibility in doped colloidal suspensions, Physical Review B 75, 224203 (2007a).
  * Ghofraniha  _et al._ [2007b] N. Ghofraniha, C. Conti, R. Di Leonardo, B. Ruzicka, and G. Ruocco, Ageing of the nonlinear optical susceptibility in soft matter, Journal of Physics: Condensed Matter 19, 205129 (2007b).
  * Mézard  _et al._ [1987] M. Mézard, G. Parisi, and M. A. Virasoro,  _Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications_ , Vol. 9 (World Scientific Publishing Company, 1987).
  * Lundgren  _et al._ [1983] L. Lundgren, P. Svedlindh, P. Nordblad, and O. Beckman, Dynamics of the relaxation-time spectrum in a cumn spin-glass, Physical review letters 51, 911 (1983).
  * Cugliandolo and Kurchan [1993] L. F. Cugliandolo and J. Kurchan, Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model, Physical Review Letters 71, 173 (1993).
  * Vincent  _et al._ [1997] E. Vincent, J. Hammann, M. Ocio, J.-P. Bouchaud, and L. F. Cugliandolo, Slow dynamics and aging in spin glasses, in  _Complex Behaviour of Glassy Systems_ (Springer, 1997) pp. 184–219.
  * Debenedetti and Stillinger [2001] P. G. Debenedetti and F. H. Stillinger, Supercooled liquids and the glass transition, Nature 410, 259 (2001).
  * Vaknin  _et al._ [2001] A. Vaknin, M. Pollak, and Z. Ovadyahu, Memory and aging in an electron glass, SPRINGER PROCEEDINGS IN PHYSICS 87, 995 (2001).
  * Mehta [2004] M. L. Mehta,  _Random matrices_ (Elsevier, 2004).
  * Glauber [1963] R. J. Glauber, Time-dependent statistics of the ising model, Journal of mathematical physics 4, 294 (1963).
  * Sherrington and Kirkpatrick [1975] D. Sherrington and S. Kirkpatrick, Solvable model of a spin-glass, Physical review letters 35, 1792 (1975).
  * Palmer [1982] R. Palmer, Broken ergodicity, Advances in Physics 31, 669 (1982).
  * Kinzel [1986] W. Kinzel, Remanent magnetization of the infinite-range ising spin glass, Physical Review B 33, 5086 (1986).
  * Henkel and Kinzel [1987] R. Henkel and W. Kinzel, Metastable states of the sk model of spin glasses, Journal of Physics A: Mathematical and General 20, L727 (1987).
  * Sibani and Hoffmann [1989] P. Sibani and K. H. Hoffmann, Hierarchical models for aging and relaxation of spin glasses, Physical review letters 63, 2853 (1989).
  * Kohring and Schreckenberg [1991] G. Kohring and M. Schreckenberg, Numerical studies of the spin-flip dynamics in the sk-model, Journal de Physique I 1, 1087 (1991).
  * Bouchaud [1992] J.-P. Bouchaud, Weak ergodicity breaking and aging in disordered systems, Journal de Physique I 2, 1705 (1992).
  * Parisi and Ritort [1993] G. Parisi and F. Ritort, The remanent magnetization in spin-glass models, Journal de Physique I 3, 969 (1993).
  * Cugliandolo and Kurchan [1994] L. F. Cugliandolo and J. Kurchan, On the out-of-equilibrium relaxation of the sherrington-kirkpatrick model, Journal of Physics A: Mathematical and General 27, 5749 (1994).
  * Cugliandolo  _et al._ [1994] L. Cugliandolo, J. Kurchan, and F. Ritort, Evidence of aging in spin-glass mean-field models, Physical Review B 49, 6331 (1994).
  * Scharnagl  _et al._ [1995] A. Scharnagl, M. Opper, and W. Kinzel, On the relaxation of infinite-range spin glasses, Journal of Physics A: Mathematical and General 28, 5721 (1995).
  * Yoshino  _et al._ [1997] H. Yoshino, K. Hukushima, and H. Takayama, Relaxational modes and aging in the glauber dynamics of the sherrington-kirkpatrick model, Progress of Theoretical Physics Supplement 126, 107 (1997).
  * Eissfeller and Opper [1994] H. Eissfeller and M. Opper, Mean-field monte carlo approach to the sherrington-kirkpatrick model with asymmetric couplings, Physical Review E 50, 709 (1994).
  * Bastolla and Parisi [1998] U. Bastolla and G. Parisi, Relevant elements, magnetization and dynamical properties in kauffman networks: A numerical study, Physica D: Nonlinear Phenomena 115, 203 (1998).
  * Nagar  _et al._ [2021] N. Nagar, N. Ecker, G. Loewenthal, O. Avram, D. Ben-Meir, D. Biran, E. Ron, and T. Pupko, Harnessing machine learning to unravel protein degradation in escherichia coli, Msystems 6, e01296 (2021).
  * Uhlenbeck and Ornstein [1930] G. E. Uhlenbeck and L. S. Ornstein, On the theory of the brownian motion, Physical review 36, 823 (1930).
  * Nyberg  _et al._ [2016] M. Nyberg, T. Ambjörnsson, and L. Lizana, A simple method to calculate first-passage time densities with arbitrary initial conditions, New Journal of Physics 18, 063019 (2016).