Exploring Complex Dynamical Systems via Nonconvex Optimization

With highly nonlinear differential models and a simple Euler integration scheme, we must be careful to ensure that the results of the optimizations are still reflective of the physicochemical model and not simply artifacts of numerical integration error. We minimize the effect of these errors in several ways.

First, we control the step size $\Delta\bm{X}_{t}$ via penalty terms in the loss. In 3.1 this is a natural part of the loss (Eq. 8) which encourages stable pattern formation, and we set $\lambda_{\Delta}=100$ in those experiments. For the remaining experiments we introduce a penalty only for step sizes that exceed a threshold $\delta_{Max}$:

Where we used $\lambda_{\delta}=1000$ and $\delta_{Max}=.05$ for all experiments except those in 3.3 which used $\lambda_{\delta}=1.0$ and $\delta_{Max}=.3$. The full loss then is the sum of $\mathcal{L}_{\delta}$ and the loss described for each experiment.

Additionally, for all results presented we verified that the results were unchanged if we ran forward simulations with the time step $\Delta t$ decreased by two orders of magnitude (from 1.0 to .01). This decreased time step should dramatically reduce numerical error and so the fact that our results are unchanged indicates they are unlikely to be artifactual.

When the number of chemical species to be visualized $N_{s}$ is $>$3 we first perform a principal component analysis (PCA) projection of the concentration dimensions. We retain only the largest 3 components and map these to the intensities of the red, green and blue channels, respectively, of an RGB image for visualization. In all cases presented here these three components explained $>$ 90% of the concentration variance. Note that this approach allows approximate visualization of arbitrarily complex chemical species mixtures, but precludes inclusion of a concentration color scale bar. This is however unimportant, given the fact that both the concentrations and reaction rate constants are determined simultaneously in the optimization; they could both be arbitrarily re-scaled and so we are effectively working in arbitrary concentration units.

We used an isotropic discrete approximation of the Laplacian operator [54] to approximate Eq. (5). In all cases diffusion coefficients were constrained to lie between $.05$ and $.02\times 10^{-5}\mathrm{m^{2}/s}$. The lower bound serves to prevent runaway accumulation of chemical species and spatial discontinuities in concentration. The upper bound was chosen to prevent numerical artifacts (checkerboard patterns, oscillation) given the $1\mathrm{s}$ time step and discrete Laplacian approximation. Using these units for the diffusion coefficients the physical dimension of a single spatial element of $\bm{X}$ (the ‘pixel size’) is $.01\mathrm{m}$, but given that the diffusion coefficients and concentrations are all optimized simultaneously and could be arbitrarily re-scaled, all the units are effectively arbitrary. All spatial domains were $64\times 64$ pixels except for in 3.3 where the dimensions were $48\times 96$.

The generator model $G(\boldsymbol{z},\boldsymbol{\theta}_{G})$ is a simple architecture consisting of repeated blocks of stride 2 convolution transpose, batch normalization [55] and a $tanh$ activation. The final generated spatial domain is therefore $2^{N_{blocks}}$ in both width and height, and the number of blocks varies within the presented results accordingly. The output is also filtered with a fixed 2D gaussian kernel ($\sigma=1.0$ pixels), to avoid introducing spatial discontinuities and therefore numerical error in diffusion simulations and diffusion entropy rate calculations. Finally, we exclude batch norm from the last layer (as in DCGAN [38]) and instead add only an optimizable scaling that is constrained to prevent negative concentrations and limit the maximum generated concentration to 10.0 (in arbitrary units).

The encoder models used in 3.2.1 use a similar structure, with stride 2 convolutions followed by a $tanh$ activation and batch normalization.

Experimentation was not extensive and constrained by hardware limitations, but in general we saw no obvious improvement in convergence properties or final optimized values from using ReLU activations or higher capacity generator or encoder architectures.

The dynamics shown in Figure 1 are derived from a ‘coupled Gray-Scott’ reaction system, which consists of two standard Gray-Scott reaction systems [56] linked by an additional reaction which allows interconversion of the two autocatalysts.

The reaction networks used here are “dense” in the sense that they contain all possible reactions matching a particular prototype. That is, if a reaction network contained species $\mathcal{S}=\{A,B,C,D\}$ then for the first reaction in (3) we would include:

While for the second reaction $A\rightleftarrows B$ we would include every possible reversible unimolecular interconversion e.g. $A\rightleftarrows B$, $A\rightleftarrows C$, $A\rightleftarrows D$, $B\rightleftarrows C$, and so on. More formally, we would generate reaction sets matching each of the reaction prototypes given in (3):

Every reaction in these sets is governed by independent forward and reverse reaction rates which are free parameters determined during the optimization. This gives e.g. a total of 40 reactions for a 5-species dense CRN (including forward and reverse), and 24 reactions for a 4 species dense CRN.

To model the dynamics of these chemical reactions we assume standard mass-action kinetics. Using the matrix representations commonly used in chemical reaction network theory, the reaction rate calculations $\Delta\bm{X}_{t_{Rxn}}$ are matrix multiplications [16][57], which are efficiently calculated on the GPUs used for both forward simulation and optimization.

We chose the reaction network described above both for the variety of reactions which comprise it, as well as it’s potential for diverse behavior. This gives the optimization process a broad space of possible reaction sets and dynamics to explore, allowing for the possibility of interesting behavior without pre-specifying it.

The reactions include many if not most stereotypical uni- bi- and tri-molecular reactions. Given sufficient timescale, more complex reaction mechanisms can be approximated with these more elementary reactions. Of more specific interest they include, as a subset, the reactions from well-studied reaction-diffusion models exhibiting complex behavior such as the Gray-Scott reaction system [56] or the Brussellator [4]. It is therefore not surprising that chemical reaction network theory tells us, on the basis of the topology and stoichiometry of these reactions, that complex dynamics are at least possible [16].

The optimization process effectively chooses a specific reaction network from the many possible networks the model described above can represent. We can visualize this reaction network as a graph, with arrows pointing from reactants to products, and with the style of arrow indicating the ‘strength’ of that reaction: Darker, larger arrows indicate larger reaction rate constants, and the graph is laid out such that strongly reacting species should be closer together (shorter arrows).

This allows us to compare the optimized reaction network from, for example, the dissipation-maximizing loss used in 3.2 (Figure 6a), to that from the replication-like-dynamics-inducing loss used in 3.3 (Figure 6b). It is worth noting that these are just the reaction rate constants, and the net flux through these reaction pathways depends on the state $\bm{X}$ and its time progression.

As is commonly observed in neural network optimization, we found that initialization was important for convergence. Specifically, in our case models initialized with uniform random but low reaction rates ($k_{fi},k_{ri}<1\times 10^{-3}\hskip 5.69046pt\forall\hskip 5.69046pti)$ and similarly low drive flow rates and concentrations failed to converge. Models initialized to more highly dissipative regimes - with higher reaction rates and higher drive flow rates - converged more reliably. Exhaustive exploration of initialization dependence is prohibitively computationally intensive, so we instead mimicked the kinetic parameters used in [1], generalizing it to our dense CRNs: One autocatalytic reaction per chemical species was initialized with a rate constant of 1. All other reaction rates were set to $1\times 10^{-3}+\epsilon_{r}$ with $\epsilon_{r}\sim\mathrm{Unif}\left(-1\times 10^{-4},1\times 10^{-4}\right)$ as this was found to be low enough to prevent large $\Delta\bm{X}$ and the associated numerical inaccuracies and/or runaway reaction rates. The small uniform random noise $\epsilon_{r}$ breaks symmetry in the reaction rates.

Feed concentrations $x_{i}$ for the flow-reactor drive were initialized to $\frac{1}{4^{N_{a_{i}}}}$ where $N_{a_{i}}$ is the number of autocatalytic reactions producing species $i$. This scaling was found to produce a highly dissipative yet somewhat kinetically and numerically stable initial state.

Diffusion coefficients were initialized from $\mathrm{Unif}\left(.05,.2\right)\times 10^{-5}\mathrm{m^{2}/s}$, uniformly within the numerically stable regime.

We used the ADAM optimizer [59], with learning rates and moving averages $\beta_{1}$ as given in Table 1 ($\beta_{2}$ of .999 was used in all experiments). With these highly nonlinear dynamical models we found that, especially for optimizations with longer timescales, gradient clipping was essential for stability of convergence. We used gradient norm clipping at 0.5 as implemented in Keras [60]. Learning rates were automatically decreased and optimization was automatically terminated via Keras’ callbacks with patience parameters as given in Table 1. In some experiments we used RMSD instead of the L2 norm given in Eq. 7 to simplify changing the size of the spatial domain without altering the magnitude of the loss, and found this gave similar results to the L2 norm as expected.

Calculating gradients with these models requires backpropagation through the entire time series $\mathcal{X}$, introducing significant memory overhead. For this reason, and because in most cases we are optimizing for a single $\bm{X_{0}}$, we set the batch size to 1 except for in Sect. 3.2.1, where we used 2. For experiments where dissipation losses were included, $\lambda_{\sigma}$ was set to 1.0 except for in 3.2.1 where it was 0.04. During optimization the reaction rate constants and feed concentrations were constrained to be non-negative, and the flow rate was constrained to between .01 and 1.0. Diffusion coefficients were constrained to within their initialization range (A.3).

For the experiments shown in 3.3 we must optimize over timescales long enough for the full replication-like dynamics to occur. Optimization from a random initialization above $\sim 250$ time points (equivalent to 250 seconds) was highly unstable, in large part due to an inability to balance vanishing and exploding gradients in the highly nonlinear CRN. Convergence was finally achieved via an “incremental” optimization approach, where we first optimize at shorter timescales and then use the parameters $\boldsymbol{\theta^{*}}$ from shorter timescales to initialize optimization at longer timescales. Additionally, at the longest timescales the variance target $\beta^{*}$ seemed to be at odds with the replication fidelity loss (Eq. 13), and decreasing this value gave better parent-daughter correlations at convergence. It is worth noting that the large number of iterations required resulted in long optimization; the final result presented here required $>600\mathrm{k}$ weight updates in total for a wall clock time of more than 26 days. The hyperparameters used for each optimization that produced Figure 5 are given in Table 2. The ADAM optimizer $\beta_{1}$ was set to .95 for all optimizations.

We performed additional analyses to test the necessity and sufficiency of the optimization of each component of the models. These computational “ablation experiments” correspond to several conditions:

• •

  “Optimized” - All components optimized.

• •

  “Random $\bm{X_{0}}$” - Optimized initial conditions $\bm{X_{0}}$ replaced with normally distributed random noise with the same mean and variance as the optimized initial conditions (truncated to prevent negative concentrations).

• •

  “Initial $\boldsymbol{\theta}_{\Psi}$” - Optimized reaction network parameters and diffusion coefficients replaced with their values from the start of optimization (as described in D.1).

• •

  “Random $\boldsymbol{\theta}_{\Psi}$” - Optimized reaction network parameters and diffusion coefficients replaced with normally distributed random noise with the same mean and variance as the optimized parameters (truncated to prevent negative reaction rates or diffusion coefficients outside the numerically stable range).

The results of these experiments are shown in Figure 2c and Figure 3c with some additional shown here and described below. These results confirm that all optimized components are required for persistent correlation with the target and for persistently high dissipation rates.

In Figure 7a we show the control experiments as in the main text but for the second example (from Figure 2b). We then repeat the same analyses for both examples but with the time step now two orders of magnitude shorter, to further reinforce that the stability and fidelity is not a result of the optimization exploiting numerical error (Figure 7b-c).

In Figure 8a we show the same result as the main text but for the second example (from Figure 3b). We also run both examples at a longer timescale to demonstrate that the optimized states are stably dissipative and produce entropy at consistently higher rates than those with non-optimized components (Figure 8b-c). Note that in this case the models with mean and standard deviation-matched random kinetic parameters (“Random $\boldsymbol{\theta}_{\Psi}$”, shown in red) are highly unstable and display significant numerical integration error (as evidenced by their behavior changing with a smaller time step). This further reinforces both the necessity and effectiveness of the numerical error avoidance methods described in Section A.1.