Generic optimization by fast chaotic exploration and slow feedback fixation

Living systems adapt to various environmental conditions by changing their internal states through processes such as gene expression and epigenetic modiﬁcation. In this paper, we propose a generic mechanism for optimization that combines fast oscillatory dynamics with a slower feedback ﬁxation process. Through extensive model simulations, we demonstrate that the fast chaotic dynamics serve as a global search for optimal states, which are then ﬁxed by the slower dynamics. This mechanism becomes more effective as the number of elements is increased. We also discuss the potential relevance of this optimization mechanism to problems in artiﬁcial neural networks.

Living systems adapt to various environmental conditions by changing their internal states through processes such as gene expression and epigenetic modification.In this paper, we propose a generic mechanism for optimization that combines fast oscillatory dynamics with a slower feedback fixation process.Through extensive model simulations, we demonstrate that the fast chaotic dynamics serve as a global search for optimal states, which are then fixed by the slower dynamics.This mechanism becomes more effective as the number of elements is increased.We also discuss the potential relevance of this optimization mechanism to problems in artificial neural networks.DOI: 10.1103/PhysRevResearch.5.023017 Biological systems can generally adapt to various environmental conditions by adjusting their internal states in response to changing environmental conditions in order to optimize survival.The most well-known and studied mechanism for adaptation is the alteration of gene expression patterns through signal transduction networks in cells, which respond to external signals.These networks have evolved over generations to adapt to specific environmental conditions.
Even though such a signal transduction mechanism has been thoroughly investigated [1], it is not fully understood whether these networks can explain all forms of cellular adaptation [2][3][4][5].Braun's research on yeasts, for example, has demonstrated that they can spontaneously adapt to different conditions, including de novo conditions that their ancestors have not experienced [2].These adaptations have been observed in artificially embedded networks, even in the absence of environmental information [3].This suggests the existence of potentially generic mechanisms for adaptation that have yet to be uncovered.In this paper, we aim to explore such mechanisms, drawing inspiration from the biological process of adaptation.If we are able to uncover a generic mechanism, it can be abstracted and applied to search for an optimized state in complex systems with many degrees of freedom.
We begin by reviewing the background of theoretical studies on the search for a generic mechanism for adaptation.One such proposed mechanism is the attractor selection, which involves the selection of a state with a higher growth rate by taking advantage of noise and growth dilution [3,6,7].While this mechanism was widely applicable [8,9] still not clear whether attractors exist in gene expression dynamics that are fitted to the environment.To address this gap in understanding, an introduction of an abstract epigenetic modification process was proposed as a means of generating stable states that differ from the original gene expression dynamics [7].While this process has the potential to enhance the applicability of the mechanism, an evolutionary process to optimize the network itself, along with the introduction of an appropriate level of noise, is needed [7,10].
Here, epigenetic modifications, such as DNA methylation or histone modification, are biomolecular mechanisms [11,12] that play an important role in stabilizing cell differentiation [13][14][15][16].Although their detailed mechanisms vary among cells [12,16,17], the modification levels generally change the feasibility of gene expression and progress more slowly than expression levels [18].
By abstracting these observations, we proposed a theory, in which the interplay between oscillatory gene expression dynamics and slower epigenetic modification generates and stabilizes novel cell types.This theory offers a potential explanation for robust cell differentiation and reprogramming [19,20], supported in part by experimental reports of temporal gene expression oscillation in pluripotent cells [21][22][23] and the oscillation of epigenetic modification levels [24].
Inspired by these studies on cell differentiation based on oscillatory dynamics [19,20], we propose a generic adaptation mechanism.Instead of making detailed connections with cell biology, this mechanism is abstract and general, to search for an optimal state that gives a requested output pattern under a given condition.The mechanism combines chaotic oscillatory dynamics with slower feedback fixation inspired by an epigenetic modification process.This interplay allows for the adaptation of gene expression in response to changing conditions.
In this paper, we explore whether this interplay can be used to achieve optimization in a general sense, rather than making detailed connections with cell biology.Possible states are explored through oscillatory dynamics and, when a desired state is approached, the slower fixation works efficiently to stabilize it.We present a simple model of regulatory networks with oscillatory dynamics and conduct extensive simulations to test the ability of the model, in order to optimize its state under a variety of external conditions, if the original dynamics show sufficiently complex (chaotic) oscillation.Our results show that the model can achieve generic optimization as long as the degrees of freedom are large enough to allow for complex dynamics, without the need for evolutionary selection of specific networks, in contrast to the previous studies that adopt noiseinduced transitions among fixed-point attractors [3,6,7].We also identify the conditions necessary for successful optimization and demonstrate that the fraction of networks satisfying these condition increases with the number of units.In addition, we explore the generality and potential applications of this optimization mechanism in the context of machine learning, specifically in neural and artificial networks.
We consider a model consisting of N units with a regulatory network and slower modification.Each unit is represented by two variables: the fast variable x i and the slow variable θ i (i = 1, 2, . . ., N).As a cell model, x i represents the expression level of the ith gene and θ i corresponds to the epigenetic modification level.In the present model, x i serves as the state variable for the unit and θ i acts as a slower modification variable that serves as a threshold for the activation of the ith unit.The units activate or suppress each other according to the regulatory matrix J i j .If J i j is positive (negative), the jth unit activates (suppresses) the ith unit [25][26][27] as given by where F (z) is a monotonic function exhibiting an on-off switch, as F (z) = tanh(βz).We set β = 40, so F (z) is close to a step function.The value of x i = 1 or −1 represent full or nonexpression of the ith unit.In Eq. ( 1), −θ i acts as a threshold for the expression of the ith unit.In other words, the slower modification level θ i determines the feasibility of the ith unit expression.As θ i increases (decreases), less (more) input from other units is required for expression.For slower modification dynamics, we adopt the simplest form of reinforcement [7,19,20,28,29], where v k (t ) is positive and indicates a positive feedback between fast expression dynamics and slower modification dynamics.If the ith unit is expressed, it is more feasible to be expressed following the increase in θ i .This positive feedback is based on the previous experimental results [17,[30][31][32].
Here, v k (t ) is the timescale of the slower modification process depending on the fast x state.It is always set to be smaller than unity, meaning that the change in θ i is slower than that in x i .Within this range, v k (t ) is increased when the state is more optimal for the kth environmental condition.The degree of optimality is given by the expression pattern of output units x m (m = 1, 2, . . ., M < N ).By introducing X k m as the target (desired) expression pattern under the kth environment, the distance between x m and X k m is calculated by Then, as the distance approaches 0, the state is optimized.Now, v k (t ) is determined as where v max = 10 −1 and b = 4, unless otherwise noted [33].Note that v k (t ) takes the maximum value v max when x m is equal to X k m , in which case the state is completely optimized to the kth environment.
We herein adopt M = 5.Then, the total number of possible M-bit target patterns with −1 or 1, described as However, the present model notably has a symmetry x ↔ −x.Considering this symmetry as X k and −X k , there are 2 M /2 = 16 independent M-bit patterns {X k }.Each environmental condition k = 1, 2, . . ., 16 has the corresponding target expression pattern X k ; the optimization to each of these M-bit patterns is examined.
We adopt random regulatory networks, whose elements {J i j } are randomly assigned either ±1, 0, with equal probability.Each of optimization trials is started from the θ i = 0 state and a randomly chosen x i .Each of the targets X k (k = 1, 2, . . ., 16) is assigned for each trial, and we run the dynamics of Eqs. ( 1)-( 4) until they reach the final stationary state.If the final x m is equal (or sufficiently close) to X k m , the model can optimize its state to the kth external condition.
Figure 1 shows the time series of x i , θ i , and v k in a successful optimization to a certain external condition.In Fig. 1(a), the fast expression dynamics x i first converge to an irregularly oscillating state.Then, the oscillatory dynamics are fixed by the slower modification θ i when x m approaches X k m and the timescale v k (t ) is increased [Fig.1(b)].Finally, throughout these transient dynamics, the state is optimized to the desired target state when x m reaches X k m .The value of v(t ) varies in time at first [Fig.1(c)], reflecting the transient dynamics of x i , until a desired state is selected.Figure 1(d) shows the optimization dynamics in x space using the principal component analysis (PCA) obtained from oscillatory dynamics with θ i = 0. Next, we investigate the optimization capacity of the model.The number of environments to which optimization is achieved among 16 environmental conditions determines the capacity. Figure 2 shows the distribution of the capacity over 500 random networks with N = 100.Here, the criterion for optimization to the kth environment is that the trial in the kth environment finishes with |x m − X k m | 2 < 10 −2 .We define the capacity as the fraction of external conditions that can be optimized at least once out of three trials [34].For most of the random regulatory networks {J i j }, the model can optimize its state to more than half of the 16 conditions.
Next, we examine the optimization to multiple conditions.Figure 3(a) shows the dynamics in x for three different environmental conditions starting from the same initial condition, using the PCA obtained from oscillatory dynamics with θ i =0.In Fig. 3(a), the gray curve corresponds to the trajectory with θ i = 0, which shows chaotic dynamics.With the slower θ dynamics of Eq. ( 2), the optimal state is reached and fixed after transient (chaotic) oscillation, depending on each target condition k.In Fig. 3(b), we study how {x i (t )} with θ i = 0 and each of target patterns {X k } comes closer by introducing the inner product of x and X k given by (1/M ) M m x m X k m to characterize the distance between x and the kth target pattern.As shown in Fig. 3(b), the time series with θ i = 0 globally explores the phase space and approaches the target optimal patterns (k = 1, 12, 14 in this example); however, it cannot approach target pattern (k = 9) where the inner product remains around 0. The distribution for inner products is extended globally over [−1, 1] for the former case, but is centered around zero for the failure case [Fig.3(c)].
We next focus on how the capacity depends on x dynamics with θ i = 0. Figure 4(i  these trajectories, x dynamics with small capacity travel small portions of phase space, whereas those with large capacity travel large portions of phase space [Fig.4(ii)].In Fig. 4(i), the former has a limit cycle attractor [Fig.4(i)(a)], while the latter has a chaotic attractor with two positive Lyapunov exponents [Fig.4(i)(b)].
To examine if the global traveling of the orbit at θ i = 0 is relevant to optimization, we computed the globalness of trajectory against the target patterns {X k m }, defined as ( In Fig. 5(a), we plot the capacity against the globalness by sampling with a 0.1 bin size and averaging the capacity for each bin.As shown in Fig. 5(a), the averaged capacity monotonically increases with the globalness of the trajectory.
As shown in Figs. 3 and 4, such global traveling is supported by chaotic dynamics.We computed the Lyapunov spectra of dynamics with fixed θ i = 0 and examined how they correlate with the capacity of the system.Figure 5 shows the correlation in the capacity against the number of positive Lyapunov exponents, sampled over 500 randomly chosen networks for N = 60.The number of positive Lyapunov exponents indicates the number of directions in which tiny perturbations can expand.Figure 5(b) suggests that the capacity increases with it.
As the fraction of the network that exhibits oscillatory and (higher-dimensional) chaotic dynamics increases with the system size N, the capacity is expected to increase as well [35].In Fig. 6(a), we plot the average capacity of networks with and without oscillatory dynamics for θ i = 0. Note that for N > 70, it becomes difficult to sample networks without oscillatory dynamics sufficiently.As N increases, more networks can be optimized to perform well under almost all (=2 M /2) environmental conditions.However, the present definition of capacity, which requires at least one successful optimization for three trials, may be too lenient.We also calculate the success rate of optimization across all trials.In the computation of these N dependencies, we choose b = 6, following the result of b dependency of success rate in N = 100 as shown in Fig. 6(c).In Eq. ( 4), b is a parameter that gives resolution to examine the distance from the target state.If b is close to 0, there is no distinction between the optimal and nonoptimal states.If b is larger, fixation works only when the state is very close to the target state so that the optimization have been completed within the simulation steps (here 10 3 ).Hence, there exists optimal b (∼6), as shown in Fig. 6(c).
In this paper, we present an optimization mechanism based on fast oscillatory gene expression dynamics coupled with a slower epigenetic fixation process.The chaotic oscillatory dynamics are relevant to the search for optimal states depending on the input, and once optimal states are approached, a slower modification process fixes those states.As long as the search using oscillatory (chaotic) dynamics adequately covers the state space, this optimization mechanism works efficiently.The degree of chaos or the region of the phase space traveled by orbits is correlated with the capacity of environments to which the cell can optimize.As the number of units (degrees of freedom) increases, the fraction of networks allowing for such dynamics also increases, supporting the generality of the proposed mechanism.
We have previously reported that the interplay between gene expression oscillation and slow epigenetic feedback allows for robust cell differentiation, which is necessary for multicellular organisms [19,20].The demonstration that chaotic oscillatory dynamics and slower fixation can allow for optimization in multiple environments may provide insight into a coherent understanding of both multicellular differentiation and unicellular adaptation.Of note, the use of chaotic exploration for optimization here is distinguishable from the earlier adaptation model in Ref. [7], where the switching among fixed-point attractors by noise is adopted.In that case, the capacity of adaptation to a variety of conditions is rather limited even after the evolution of networks under given fitness.In contrast, the exploration of chaotic dynamics here allows for optimization (adaptation) to almost all conditions, without tuning the network, as long as the system is sufficiently high dimensional.
Our optimization mechanism can be depicted as Waddington's epigenetic landscape, which is often used by developmental biologists [36].The initial state in a shallow valley, in the epigenetic landscape, travels over a large portion of phase space, similar to the chaotic dynamics in our model.As the slow variables change, deep valleys are generated to which the state is attracted, leading to optimization or adaptation.
The scheme proposed in this paper does not require attractors that achieves optimization in advance or the evolutionary optimization of networks.In this sense, it would also support the generic and spontaneous adaptation of cells to unforeseen environmental conditions.While there is some support for oscillatory expression dynamics in stem cells in multicellular organisms [21][22][23]37,38] and theoretical verifications [39][40][41], there is currently no direct evidence for oscillatory expression dynamics in unicellular organisms.Expression dynamics are typically noisy and it is challenging to extract oscillatory components experimentally.Notably, our mechanism works robustly under strong stochasticity.
Our model adopts a simple setup for oscillatory dynamics with on-off type dynamics and a slower fixation process.Onoff dynamics are widely used in biological and artificial neural networks.The present scheme uses an autonomous search for the desired state through chaotic dynamics and slower fixation can be generally applied to learning or optimization processes.Compared to random sampling used in simulated annealing [42], chaotic dynamics do not need to sample the entire space, making the search more efficient [43][44][45].We also note that the relevance of chaos or chaotic itinerancy to neural information processing has been discussed [46][47][48][49], whereas θ i in our model can be regarded as inputs [50].In contrast to the Hebbian learning, which requires changes to J i j (i.e., N × N elements), our scheme requires changes only to θ i (i.e., N elements), which would be useful for low-rank changes in reservoir computation or echo-state networks [51,52].
, it wasPublished by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

FIG. 1 .
FIG. 1. Optimization process of our model with N = 100.Time series of (a-i), (a-ii) x i (i = 1, 2, . . ., N ), (b) θ i , and (c) v k .Starting from the initial condition with random x i and θ i = 0, the state converges to an oscillatory state.With gradually developed θ i , the state reaches the desired expression pattern X k m .That is, v k takes the maximum value v max (=10 −1 ).(a-ii) Time series of x i only for i = 1, 2, . . ., 5 for t ∼ 10-20.(d) Optimization dynamics in x space.We adopt PCA obtained from oscillatory dynamics with θ i = 0.The state starting from a random initial condition (red point) reaches the target gene expression pattern (blue point) throughout transient oscillatory dynamics.
FIG.2.Distribution of capacity computed from 500 random network models with N = 100.

FIG. 3 .
FIG. 3. (a) Optimization process against three different condi-(targets) plotted in x space using the PCA space for oscillatory dynamics with θ i = 0 (gray trajectory).The colored trajectories show optimization in three different environments (k = 1, 12, 14) starting from identical initial condition (x).Each optimization is completed at the colored circle.(b) Time series (left) and histogram (right) of the overlap of {x m (t )} with θ i = 0 and target patterns (1/M ) M m x m X k m .Black time series and histogram represent the case of nonachieved environments (k = 9).The trajectory with fixed θ = 0 is not influenced by the target and maintains chaotic dynamics.
FIG. 4. Comparison of the dynamics with fixed θ i = 0 for (a) small and (b) large capacity with N = 60.(i) Dynamics plotted in the PCA space of x. (ii) Time series of (1/M ) M m x m X k m and the overlap of {x i (t )} with θ i = 0 and target patterns {X k }.Left: Capacity = 3, globalness = 0.16, and no positive Lyapunov exponents.Right: Capacity = 13, globalness = 0.35, and two positive Lyapunov exponents, 0.51 and 0.16.

Figure 6 (
FIG.5.Capacity vs the characteristics of the dynamics for fixed θ i = 0. (a) Capacity as a function of the globalness of trajectories, as defined in Eq. (5).(b) Capacity as a function of the number of positive Lyapunov exponents, N = 60.