Molecular Engineering of Robustness and Resilience in Enzymatic Reaction Networks.

Living systems rely on complex networks of chemical reactions to control the concentrations of molecules in space and time. Despite the enormous complexity in biological networks, it is possible to identify network motifs that lead to functional outputs such as bistability or oscillations. One of the greatest challenges in chemistry is the creation of such functionality from chemical reactions. A key limitation is our lack of understanding of how molecular structure impacts on the dynamics of chemical reaction networks, preventing the design of networks that are robust (i.e., function in a large parameter space) and resilient (i.e., reach their out-of-equilibrium function rapidly). Here we demonstrate that reaction rates of individual reactions in the network can control the dynamics by which the system reaches limit cycle oscillations, thereby gaining information on the key parameters that govern the dynamics of these networks. We envision that these principles will be incorporated into the design of network motifs, enabling chemists to develop "molecular software" to create functional behavior in chemical systems.

Introduction

Vast metabolic and genetic networks of chemical reactions allow living cells to sense their environment, react to stimuli, and use nutrients for cell growth and division.[1] Although these networks are daunting in complexity, recurring patterns, so-called network motifs, have been identified that create functional behavior on a smaller scale.[2,3] Simple motifs with a few positive and negative feedback loops create functionality such as bistable switching, adaptation, and oscillations.[4−6] The key challenge for chemistry is to translate the design principles of living systems into robustly engineered artificial systems.[7−9] Chemical reaction networks organized into different motifs give rise to rich dynamic behavior, but programming their precise output has proven very delicate. Early work has resulted in a number of exciting examples ranging from functional out-of-equilibrium systems that can perform logic operations[10,11] to dissipative self-assembling structures creating new forms of smart materials.[12−17] However, it is clear that we do not fully understand how to engineer robustness and resilience in molecular reaction networks.[18−21]

Strategies to obtain robustness and resilience usually rely on modifying or increasing the networks’ connectivity,[22−27] but this approach fails to take into account the nature of chemical reactions within the networks. A bottom-up construction of complex molecular systems offers a novel route to directly probe the influence of molecular reactivity on the dynamics of reaction networks.[28−30] We recently reported a rationally designed, fully characterized enzymatic reaction network showing limit cycle oscillations (Figure a).[31] This two-node oscillator integrates the autocatalytic production of the enzyme trypsin with a tunable and delayed negative feedback induced by trypsin-activated small molecules. Unlike most biological systems, all rates in our network are known and can be individually addressed.[32] Here, we synthesized a small library of pro-inhibitor molecules (Figure a) to finely tune rate constants for nodes controlling the activation and the termination of the negative feedback in topologically equivalent CRNs.

Figure 1

Enzymatic reaction network with modified pro-inhibitors leading to a library of CRNs. (a) Topology of the enzymatic reaction network composed of enzymes trypsinogen (Tg), trypsin (Tr), and aminopeptidase (Ap) and a synthetic pro-inhibitor. Modifications to the pro-inhibitor were achieved by conventional synthetic procedures (see Supporting Information S1). Substituents acetyl (Ac), methoxy ethelene glycol (MEG), and acetyl aspartic acid (Ac-Asp) and amino-methyl (Me), ethyl (Et), or propyl (Pr)-benzenesulfonyl fluoride were used for R1 and R2, respectively. (b) Illustration of CRN matrix composed of nine different pro-inhibitors combined from R1 × R2. The rate constants for activation, kact, and inhibition, kinh, were determined in isolated reactions (see Supporting Information S1).

We show how the steepness of the response of our negative feedback can be controlled by the substituents in the pro-inhibitor molecule and how this approach leads to more robust networks that reach their stable limit cycle more quickly. It is exactly this understanding of the dynamics of the networks that allows us to identify the chemical principles that govern how complex systems reach their out-of-equilibrium function.

Results and Discussion

Our network combines a positive and a delayed negative feedback loop. In the reaction network, trypsin (Tr) catalyzes its own formation from the precursor trypsinogen (Tg). Opposed to this positive feedback, Tr is inhibited by the negative feedback that is composed of three sequential steps (Figure a). In the activation step, Tr converts a pro-inhibitor into an intermediate inhibitor (Int-I), which consists of a glutamine (Gln) residue attached to a potent inhibitor for Tr. Another enzyme, aminopeptidase N (Ap), controls the release of the inhibitor moiety by cleaving Gln in the delay step. In the final step, Tr recognition of the active inhibitor closes the negative feedback loop.

We tune the rates in the activation and inhibition steps by modifying the pro-inhibitor structure. The combination of three substituents on R1 and R2 sites results in nine pro-inhibitors with different reaction rate constants for either the activation or final inhibition step of the negative feedback loop (Figure b). Substituents acetyl (Ac), methoxy ethylene glycol (MEG), and acetyl aspartic acid (Ac-Asp) on position R1 influence the rate constant of activation (kact) of the negative feedback loop by changing the affinity of the pro-inhibitor toward Tr.[32,33]

On the other end of the molecule, changing the length of the alkyl chain (methyl, ethyl, and propyl) of the inhibitor in position R2 allows us to tune the rate constant of inhibition (kinh) by changing the actual fit of the inhibitor in the active pocket of Tr.[34,35] Overall, we obtained a small family of networks (chemical reaction networks (CRNs) 1–9) with the same topology, but with different kinetics, as indicated by the different lengths of arrows in Figure b (the shorter the arrow, the stronger the interaction and the higher the reaction rate constant).

The pro-inhibitor (Pro-I) is the source of the negative feedback and is therefore responsible for carrying the reaction back to its original state. As evidence for the fine-tuning of the negative feedback, we first confirm in batch experiments that each CRN containing a different pro-inhibitor can exhibit the desired function. The response of the reaction network in batch conditions gives an initial rise in [Tr] before decaying to equilibrium in which the active inhibitor is bound irreversibly to Tr. The series of experiments in Figure a show that all networks exhibit similar behavior, but also demonstrate that the subtle changes in the pro-inhibitor molecules have changed the details of the response.

Figure 2

Network properties in equilibrium conditions. (a) Experiments with CRNs 1–9 carried out in identical batch conditions (i.e., with initial conditions [Tg]0 = 130 μM, [Pro-I]0 = 260 μM, [Tr]0 = 0.2 μM, [Ap]0 = 0.8 U mL–1). Trypsin concentration is determined by a fluorescent assay (see Methods summary). (b) The inhibition time (purple domain in (a)) is defined as t(Tr,max) – t(Tr=0.2) and plotted as a function of kinh (c) Experiments with CRN series 4, 5, 6 and 3, 6, 9 carried out in identical batch conditions in the absence of Tg ([Pro-I]0 = 200 μM, [Tr]0 = 100 μM, [Ap]0 = 0.8 U mL–1). The concentrations of methyl- to ethyl- to propylamine inhibitors (4, 5, 6) are depicted in gray, light blue, and dark blue squares, respectively. Inhibitor concentrations are determined by UV detection after being separated by HPLC (see Methods summary). (d) The response of the negative feedback determined the gradient (d[Tr]/d[Act-I]) in which [Tr] decreases from 80 μM to 20 μM (see Supporting Information S2).[36]

There are clear trends of gradual changes in peak position, area under the peak, and maximum [Tr] obtained. We analyzed the peak characteristics of individual responses in detail (see Supporting Information S2), but here we wish to highlight (in purple) how the modifications on Pro-I influence the time required to bring the maximum [Tr] from the peak back to zero. This “inhibition time” gives information on the “strength” of the negative feedback. We find that the inhibition time correlates negatively to the inhibition rate constant kinh, and this trend is consistent for the three different substituents in R1 (Figure b).

The library of CRNs also provides deeper insights into the nature of the negative feedback. Choosing two series of CRNs allows us to experimentally isolate the influences of the various rate constants in our enzymatic reaction network on the negative feedback loop. We used CRNs 4, 5, and 6 (with various substituents for R2) to investigate the impact of changes in the inhibition rate constant (kinh) and CRNs 3, 6, and 9 (with various substituents for R1) for changes in the activation rate constant (kact). Figure c shows the various responses of the negative feedback all starting with initial trypsin concentrations [Tr]0 = 100 μM. In the absence of Tg, the negative feedback initiates immediately, as can be seen in the decay in [Tr] and the simultaneous increase observed in [Act-I]. Figure c shows Tr is inhibited faster in the series 6 > 5 > 4 (requiring less time to reach [Tr] = 0). Furthermore, the Act-I production during the reaction is significantly slower in the same series. In contrast, changes in the kact (CRNs 3, 6, and 9) do not influence the time to fully inhibit Tr nor the amount of Act-I that is produced in the reaction.

The effect of kinh and kact on the negative feedback reaction is studied in more detail by following how [Tr] changes as a function of [Act-I]. To evaluate the kinetic interplay between the key enzyme Tr and the active inhibitor that is eventually produced in the solution, we determined the gradient (d[Tr]/d[Act-I], Figure d). We note that this gradient changes from a gradual (dTr/dAct-I = −0.99) to a much steeper, so-called ultrasensitive,[36] response (dTr/dAct-I = −4.41) when changing from methyl- to ethyl- to propylamine inhibitor (4, 5, 6) (e.g., higher rate constant kinh) in the network. This is important, as proper balancing of the time scales of opposing chemical reactions is necessary in order to obtain sustained oscillatory behavior under out-of-equilibrium conditions.[6] Consequently, a gradual decay in [Tr] as a function of [Act-I] is expected to lead to a negative feedback loop that would be less to counteract the autocatalytic production of Tr in a timely fashion. Importantly, Figure d demonstrates that the steepness of the gradient responds only to changes in R2 (kinh), as the CRNs with changes in R1 (Ac, MEG, or Ac-Asp groups on R1) all show the same response. Thus, the kinetics of the feedback loop are dominated by the structure of the active inhibitor.

Next, we used simulations to investigate the effect of increasing the strength of the negative feedback on the robustness of the networks under out-of-equilibrium conditions.[19,20]Figure a illustrates how we use flow (with flow rate constant kf)[37] to maintain the network out-of-equilibrium. In flow, our system exhibits sustained oscillations only in a limited parameter space. Hence, we first determined (using our previously published numerical search method)[32] the range of feed concentrations, as well as flow rates, that will lead to sustained oscillations for each of the CRNs.

Figure 3

Robustness in out-of-equilibrium enzymatic reaction networks. (a) Network motif with a continuous influx of reactants ([Tg]0, [Ap]0, [Pro-I]0). (b) Phase plot showing the predicted conditions for sustained oscillations. We identified a value of [Pro-I]0 for which the region of (kf, [Ap]0)-space (in which sustained oscillations are observed) is largest (indicated by [Pro-I]0opt). (c) Phase plots at the optimal conditions for [Pro-I]0 as a function of the parameters [Ap]0 and kf. The flow rate constant kf (in h–1) is obtained by dividing the applied flow rate (in μL h–1) over the volume of the reactor (in μL) in the simulations. (d) Normalized areas, in which the simulations predict oscillations, as a function of both kinh (differences along dashed lines) and kact (differences among colored squares).

Figure b shows the parameter window (composed from the feed concentrations [Ap]0, [Pro-I]0, and flow rate) as a “volume” of the oscillatory regime.[38] Using a computer script, we took slices of the 3D plot in the (kf, [Ap])-planes to determine the apparent optimal [Pro-I]0opt (i.e., the feed concentration of the pro-inhibitor giving the largest parameter space with sustained oscillations; see Supporting Information S3.2.2). If we repeat this procedure for all CRNs, we find that there are significant differences in the size of parameter space in which sustained oscillations can be found (Figure c). We compared the areas of the parameter space that leads to sustained oscillations for each CRN in Figure d. In contrast to the results in batch conditions, the robustness of the oscillatory regime increases to similar degrees as a function of both kinh (differences along dashed lines) and kact (differences among colored squares). Nonetheless, we find that the network is most robust with the propyl derivative of the active inhibitor.

Of particular interest is not the size of the oscillatory regime but how fast our system reaches the stable limit cycle.[39,40]Figure a illustrates in a phase portrait how the network approaches the limit cycle from various points. The number of orbits that is required to reach the limit cycle is a first indication of the resilience of the network, as the system attempts to “recover” the “perturbed” states.[41] The system’s recovery from a perturbation determines how sustained oscillations in Figure b are established after a certain period in which the higher amplitudes in the beginning of the experiments relax to their limit cycle values. Hence, this decay essentially provides information on the attractor strength of the network that “pulls” the network into oscillations.

Figure 4

Resilience in enzymatic reaction networks. (a) Illustration of a system approaching a limit cycle as the programmed behavior of our network. (b) Illustration of the fitting algorithm applied to a time trace showing oscillations in [Tr](t). The computer algorithm recognizes the sustained oscillations when at least three peaks have identical amplitudes (indicated by a purple background). If such a case exists, the script identifies the local maxima of the response (and fits them with an exponential function of the form f(t) = a0 + a1 exp(bt). (c) Comparison of attractor landscapes of CRNs 1–3 at their respective [Pro-I]0opt, based on calculated exponential decays in the relaxation period explained in (b) (see Supporting Information S4 for more details). The [Pro-I]0opt used in these landscapes are [1]0opt = 1.70 mM, [2]0opt = 1.50 mM, and [3]0opt = 0.80 mM. Each grid in the phase plots represents a simulation of 150 h with the steepness of the decay (e.g., the magnitude of the exponent) indicated by the color bar. (d–f) Experiments with CRNs 1–3 in flow conditions were carried out in a continuously stirred tank reactor with an internal volume of 118.2 μL. Six glass syringes were loaded: 4 containing [Tg]0, [Tr]0, [Ap]0, [Pro-I]0 and 2 containing a buffer solution and a fluorogenic substrate, respectively (see Supporting Information S4). The experiments are carried out with [Tg]0 = 167 μM, [Tr]0 = 0.2 μM, and [Ap]0, [Pro-I]0 and kf as reported in the graphs. See Methods summary for experimental details.

The decay toward sustained oscillations can be quantified using a mathematical model that we have adapted from our previous studies. We implemented an additional algorithm that locates the local maxima of the oscillations in the simulations. Subsequently, we fit the maxima in the relaxation period with an exponential function (f(t) = a0 + a1 exp(bt)) as depicted in Figure b. The magnitude of the exponent that results from the fitting algorithm can then be plotted as a function of [Ap]0 and flow (see Supporting Information S3.2.3), creating “attractor landscapes” in the previously found [Pro-I]0opt as shown for networks 1–3 (Figure c).

We selected the network series 1–3 to elucidate the correlation between the relaxation dynamics and the modifications on R2 in the Pro-I. The landscapes in Figure c characterize how fast the limit cycle of the networks is approached and clearly show the considerable enlargement and deepening of the regions with steep decays within the series Me, Et, Pr. The broader basin of attraction in combination with the larger exponents indicates that networks consisting of the propyl inhibitors are more resilient. We experimentally validated the trend found by numerical simulations, by studying the CRNs in flow, using feed concentrations [Ap]0 and the flow rate predicted by the simulations in Figure c. Our previous work on CRN 2(31) (Figure e) showed that this network reached sustained oscillations after approximately four oscillations. In comparison, Figure d–f shows that CRNs 1–3 all exhibit oscillatory behavior, but the transition time required to reach the desired oscillations is significantly reduced within the series. Remarkably, CRN 3, containing a propyl inhibitor, establishes sustained oscillations almost immediately (after the first oscillation). In contrast, CRN 1, containing a methyl inhibitor, requires at least seven oscillations to reach the sustained oscillations.

To verify that an enhanced resilience can be attributed to the propyl inhibitor, we performed a similar experiment using CRN 9. This control experiment shows that sustained oscillations are also reached significantly faster (requiring only two oscillations), indicating that more resilient networks can be engineered by tuning of R2 on the Pro-I (Supporting Information Figure S4.3). Notably, the trend of increasing resilience with larger kinh coincides with the same trend observed for the increase in robustness with larger kinh and can all be traced back by the changes in kinetics of the negative feedback loop as shown in Figure f.

Conclusions

We have shown how careful tuning of the reactivity of different parts of small molecules allows us to systematically engineer the responsiveness of enzymatic reaction networks. Using oscillating networks as a paradigm, we show that molecular engineering can be used to change the steepness of the response in the negative feedback loop, thus creating networks that not only show sustained oscillations over a larger parameter space, but also reach these oscillations much more rapidly. Our experimental observations are supported by a detailed computational analysis of the networks and identify the principles that govern how molecular structure impacts the dynamics of out-of-equilibrium systems. These studies pave the way for the future forward engineering of robust and resilient functional out-of-equilibrium systems using “molecular software”. Creative application of synthetic chemistry can be used to create a range of new chemical reaction networks with desired functional outputs. We believe that a focus on the strength of the individual, local interactions between components in a network will also lead to a further understanding of the functioning of biological networks.[8,24,42]

Methods

Full details of the synthesis and characterization of all compounds, kinetic studies, computational simulations, and flow experiments appear in the Supporting Information.

Batch Experiments

For the batch experiments of the full networks, various pro-inhibitors (260 μM) were mixed independently with trypsinogen (130 μM), trypsin (0.2 μM), and aminopeptidase (0.830 U/mL) in 100 mM Tris buffer, pH 7.7, containing 20 mM CaCl2. For the batch experiments with the isolated negative feedback, various pro-inhibitors (200 μM) were mixed independently with Tr (100 μM) and aminopeptidase (0.830 U/mL) in 100 mM Tris buffer, pH 7.7, containing 20 mM CaCl2. Aliquots were taken from the reaction mixture to monitor trypsin activity by a fluorogenic assay (vide infra) and inhibitor concentration by an HPLC analysis (vide infra).

Trypsin Activity Assay

Trypsin activity was measured by mixing 100 μL of the quenched reaction mixture with 3 mL of 5 μg/mL bis(Cbz-l-Arg)-rhodamine fluorogenic substrate in 50 mM Tris-HCl, pH 7.7. The increase in fluorescence intensity (λex = 450 nm, λem = 520 nm) was monitored for 40 s, and the initial, linear slope was compared to a calibration curve to find the concentration of active trypsin (see Supporting Information S2).

Determination of Inhibitor Species Concentration

A 140 μL amount of the quenched reaction mixture was filtered to remove all enzymes. The organic compounds in the filtrate were separated by analytical HPLC and were monitored in time with UV detection at 265 nm. Appropriate peaks were integrated, and a calibration curve was used to determine the concentration of inhibitor species. The calibration curve is provided in the Supporting Information S2.

Flow Experiments

Four glass syringes were loaded with trypsinogen (8 mg/mL, 338 μM in 4 mM HCl, 36 mM CaCl2), trypsin (27 μg/mL, 1.16 μM in 500 mM Tris-HCl, 20.5 mM CaCl2, pH 7.7), pro-inhibitor (5 times the desired final concentration in the CSTR, which varies, in 2 mM HCl), and aminopeptidase (10 times the desired final concentration in the CSTR, which varies, in 10 mM Tris-HCl, 10 mM MgCl2, pH 7.7) and connected with tubing to the four inlets of a 118.2 μL polydimethylsiloxane reactor. Typical reactor flow rates lie in the range of 20–35 μL h–1. Fractions of the total flow rate were 0.5 for trypsinogen, 0.2 for both trypsin and pro-inhibitor, and 0.1 for aminopeptidase. Subsequently, these fractions multiplied by the syringe concentrations determine the feed concentrations of the reactor as reported in Figure .

Two additional glass syringes were loaded with a buffer solution (50 mM Tris-HCl, pH 7.7) and a fluorogenic substrate solution (25 μM bis(Cbz-l-Arg)-rhodamine in 500 μM HCl). First, the outflow of the reactor was coupled to a mixing chamber of 10 μL, in which the content was diluted with the buffer solution. Subsequently, the diluted reaction content was connected to a dolomite T-junction chip and mixed with the fluorogenic substrate solution, to monitor trypsin activity online. Fractions of the buffer and fluorogenic solution flow rates depend on the out-flow rate of the reactor, which varies for each experiment (see Supporting Information S4 for further details). Fluorescence intensity (at λex = 470 nm, λem = 525 nm) was monitored on an Olympus IX81 inverted microscope.

Computation

At the core of all our simulations, trajectories of the individual species are simulated by numerical integration from an initial state of the system. We analyzed the key characteristics of simulated responses to identify and classify the steady states using a classification algorithm written in Matlab.[31] The overall response is considered a sustained oscillation when at least three consecutive peaks show no difference (within a defined confidence interval of 97%) to their neighbors. The time series in the algorithm are ran for 300 h.

To identify the optimal feed concentrations, [Pro-I]0opt, we calculated the probability (i.e., size of oscillatory regime in the (kf, [Ap])-planes) to maintain the desired behavior of the network (i.e., oscillations) as normalized volumes for each network.[32] For the construction of attractor plots, we wrote an algorithm that identifies and fits the local maxima in the relaxation period. The resulting exponents were summarized and plotted using a surf function in Matlab to obtain Figure c. See Supporting Information for more details.