Parasitic behavior in competing chemically fueled reaction cycles

Non-equilibrium, fuel-driven reaction cycles serve as model systems of the intricate reaction networks of life. Rich and dynamic behavior is observed when reaction cycles regulate assembly processes, such as phase separation. However, it remains unclear how the interplay between multiple reaction cycles affects the success of emergent assemblies. To tackle this question, we created a library of molecules that compete for a common fuel that transiently activates products. Often, the competition for fuel implies that a competitor decreases the lifetime of these products. However, in cases where the transient competitor product can phase-separate, such a competitor can increase the survival time of one product. Moreover, in the presence of oscillatory fueling, the same mechanism reduces variations in the product concentration while the concentration variations of the competitor product are enhanced. Like a parasite, the product benefits from the protection of the host against deactivation and increases its robustness against fuel variations at the expense of the robustness of the host. Such a parasitic behavior in multiple fuel-driven reaction cycles represents a lifelike trait, paving the way for the bottom-up design of synthetic life.

In the Supporting Information, we provide Supporting Tables with the characterization of succinate derivatives and anhydrides and with the concentrations of the anhydrides in the aqueous and oil droplet phase, respectively which is used for the construction of phase diagram. We describe our kinetic model and compare the experimental data with the theory. We tabulate the reaction rate constants as obtained from fitting routines. We also calculate the droplet composition for different competitor 2 concentrations. Furthermore, we provide additional data for the co-phase separation mechanism, periodic fueling and buffer capacity.

Supporting Tables
Supporting Table S1: Characterization of succinate derivatives (xx = detection not possible).

Theoretical kinetic model for co-phase separation and chemical reactions competing for fuel Fuel driven chemical reactions
Our system is composed of four different types of molecules: fuel, succinate derivatives, intermediate and anhydride molecules. The fuel molecule EDC is abbreviated by the letter F . We consider three types of succinate derivatives which we refer to as precursor (A P ), competitor 1 (A C 1 ) and competitor 2 (A C 2 ), respectively, and abbreviate each derivative as A i with i = P, C 1 , C 2 . The succinate derivatives used experimentally are (E/Z)-2-buten-1-ylsuccinic acid (precursor), succinic acid (competitor 1) and (E/Z)-2-hexen-1-ylsuccinic acid (competitor 2). The intermediate molecule, AF i , is composed of fuel and the respective succinate derivative. There are three corresponding anhydrides, abbreviated as B i , and referred to as product, product of competitor 1 or product of competitor 2. The anhydrides used experimentally are (E/Z)-2-buten-1-ylsuccinic anhydride (product), succinic anhydride (product of competitor 1) and (E/Z)-2-hexen-1-ylsuccinic anhydride (product of competitor 2). In our systems, the solvent is water (w). The fuel F can undergo two chemical reactions: F can get slowly hydrolysed to waste W with a rate constant k 0 (reaction 0, Eq. (1a)), or, the fuel F drives the transition from the succinate derivative A i to the intermediate AF i with a chemical flux that is proportional to the fuel concentration c F with a rate constant k f,i (reaction 1, second order chemical reaction). This intermediate molecule, AF i can spontaneously hydrolyse back to A i with a rate constant k 3,i , or irreversibly turn over to the anhydride B i (reactions 1-3, Eq. (1b)). The turn-over to the anhydride (reaction 2) occurs spontaneously as an intramolecular reaction with a rate constant k 2,i (first-order reaction). In aqueous media, this anhydride hydrolyses and thereby turns over to the initial succinate derivative (Eq. (1c)), which we shortly refer to as deactivation step. Considering a constant pH and approximately dilute conditions relative to water, deactivation follows a first order chemical reaction with a deactivation rate constant k d,i . All reactions (0)-(4) as shown in Fig.1, can be summarised by the following reaction schemes: Reaction 4:

Co-phase separation of anhydride molecules
At concentrations larger than their respective solubilities (Supporting Table S4, S5), each anhydride B i (product and product of competitor 2) can phase separate from the water-rich solvent w and form oil-like droplets that are rich in anhydride. Since fuel F , the succinate derivatives A i (precursor and competitor 2) and intermediate molecules AF i , are well soluble in water, the droplets dominantly contain the rather hydrophobic anhydrides B i with concentrations inside compared to the outside similar to oil drops in water. Here, we study the interplay of the product (P ) and the product of competitor 2 (C 2 ), which co-phase separate. To estimate the co-phase separation properties of these two anhydrides with respect to the solvent, we consider the limit of fuel excess. In this case, the system is mostly composed of the two anhydrides, solvent and fuel. Due to the hydrophilic property of the fuel, we neglect its effects on phase separation. Thus, we can determine a ternary phase diagram for the remaining three molecules (product and product of competitor 2, and solvent w). As a model for this phase diagram, we consider a ternary, incompressible Flory-Huggins free energy density, where the summation index i, j = P, C 2 , w, runs over all three types of molecules, i.e., product, product of competitor 2 and water as the solvent, respectively. Moreover, r i = ν i /ν w , where ν i is the molecular volume of component i and ν w is the molecular volume of water. Incompressibility implies that all molecular volumes are constant and that the volume fractions φ i obey, 1 − φ w = φ P + φ C 2 . These volume fractions are related to concentrations, c i = φ i /ν i . Both hetero-and homotypic interactions between molecules i and j are accounted for by the interaction parameter χ ij . Our model for the phase diagram has 5 parameters: the molecular volumes of the two anhydrides, ν P and ν C 2 , the interaction parameters between each anhydride and solvent, χ P,w and χ C 2 ,w , and the effective interaction parameter between the two anhydrides, χ P,C 2 . The ternary phase diagram is obtained from the free energy density (Eq. (2)) by a Maxwell construction, where the volume fractions in droplet phase (I), φ I,eq P and φ I,eq C 2 , coexist with the volume fractions in aqueous phase (II), φ II,eq P and φ II,eq C 2 . The volume fractions fulfil the equilibrium conditions for phase coexistence: where the exchange chemical potentialsμ i = ν i ∂f /∂φ i and the osmotic pressure Π = −f + i µ i φ i /ν i can be calculated from the free energy density (Eq. (2)), and i = P, C 2 . In the Maxwell construction for a ternary, incompressible mixture, we have two more conditions, namely the conservation laws for the total volume fraction for each anhydride, The equilibrium volume fractions in both phases for each anhydride (φ I,eq i , φ II,eq i ), together with the droplet volume V I , give us five unknowns and we have five equations determining these unknowns.

Obtaining interaction parameters from the experimental phase diagram
Using the observed molar masses (Supporting Table S1), we find for the product (P ) and the product of competitor 2 (C 2 ), m P /m w = 8.56 and m C 2 /m w = 10.11, respectively. To fit the experimental phase diagram, we used the molecular masses relative to water as initial guesses for the fractions of the molecular volumes, i.e., r i = ν i /ν w m i /m w , where m i denotes the molecular mass of molecule i. Given the experimental equilibrium concentrations of the anhydride molecules i (φ I,eq i , φ II,eq i ), in their respective binary system, i.e., product (P ) with solvent (w) and product of competitor 2 (C 2 ) with solvent (w), we can solve Eq. (3) for χ P,w , ν P , χ C 2 ,w and ν C 2 . In the main text φ II,eq is referred to as c out .
Thus, for the fraction of molecular volumes r i = ν i /ν w , we obtain r P = 6.44 and r C 2 = 8.35 which are in good agreement with the mass fraction (see previous paragraph). The obtained interaction parameters in the respective binary system are χ P,w = 1.63 and χ C 2 ,w = 1.76 (in units of k B T ). Keeping these interaction parameters and the molecular volume fractions fixed, we have only one undetermined parameter left, namely χ P,C 2 . This parameter is obtained by finding the best agreement with the binodal lines and tie line slopes. Very good agreement is obtained for the value χ P,C 2 = 0, see Fig.2. A zero χ P,C 2 parameter is consistent with the homotypic interactions and heterotypic interactions between the two anhydrides being approximately of the same magnitude; a scenario that is reasonable due to similarity of the molecular structures of the two anhydrides.

Figure 2: Ternary phase diagram.
A) The concentration of the product of competitor 2 in the aqueous phase as function of competitor 2 concentration. The value remains nearly constant around 2 mM. B) The concentration of the product in the aqueous phase as function of competitor 2 concentration. In the absence of competitor 2, the value is 27.8 mM (solubility) and it decreases with increasing competitor 2 concentration. The measurements are performed at 16 mins into the reaction cycle. C) Linear representation of the ternary phase diagram, highlighting that it spans over 3 orders of magnitude in concentration values. The circles denote the experimental data points and the solid lines represent theoretically determined binodal lines and tie lines obtained from solving Eq. (3). The dashed black lines are the experimental tie lines. In the inset, we show the equilibrium concentrations in the aqueous phase (II). The compositions of the aqueous phase and oil droplet phase differ depending on which tie line the total anhydride concentrations lie.

Kinetic model for phase separation of products and fuel-driven chemical reactions
To describe the kinetics of all the reacting molecules in the system, we introduce the concentrations of the two anhydrides (product P and product of competitor 2 C 2 ) inside the droplet phase (I) as c I B i , where i = P, C 2 and the concentrations of all other molecules outside the droplet as c j , where j = F, A i , AF i , B i denotes fuel, succinate derivates (precursor and competitor 2), intermediate molecules, and anhydrides (product and product of competitor 2), respectively. Within our model, these components can undergo chemical reactions as described in Eq. (1), and the anhydrides (product and product of competitor 2) can phase-separate. Supported by the experimental observation that droplets form very quickly on the experimentally relevant time scales of minutes, we propose a simplified model for the total concentrations within the demixed region of the phase diagram. This model relies on fast phase separation and partitioning kinetics of the anhydrides relative to their chemical reactions. Specifically, this model is valid if the inter-droplet distances do not exceed the reaction-diffusion length scale D i /k d,i for both anhydrides. Using experimental values, this reaction-diffusion length scale is in the order of a few hundreds µm, while inter-droplet distances are about a few tenths of µm, supporting the validity of this approximation. Thus, the reaction-partitioning equations for the kinetics of all molecules outside the droplets (II) and the kinetics of the total anhydride concentrations read: Under the valid assumption of fast partitioning kinetics compared to the slow chemical reactions outside the droplets, the total concentrations of the two anhydrides,c B i (t) determine their equilibrium concentrations (via Maxwell construction using Eq. (2) in Eqs. (3)) and the volume of the dense phase at each time point t: Volume dense phase: for i = P, C 2 and V is the total volume of the system which is a constant (V = 1.05 mL).
In summary, in our model which considers the limit of fast phase separation compared to chemical reactions, the chemical reactions affect the total concentrations of both anhydrides, leading to instantaneous changes of their respective equilibrium concentrations. These equilibrium concentrations varying with time modify the concentration levels in both phases (I and II) via partitioning and thereby in turn affect the chemical reactions. Due to this feedback between chemical reactions and phase separation, the total concentrations move on a specific path in the two dimensional phase diagram which is unique to the initial succinate derivative and fuel concentrations (see Fig.  3E in main text). The total concentrations of the anhydrides, i.e.,c B i (t), can also lie outside the demixed region during the reaction cycle, implying that phase separation is absent and the system is homogeneous (mixed phase). In that case, there is no partitioning and V I = 0, such that only the kinetic equations, i.e., Eqs. (4a) (4b) (4c), and Eq. (4e), are to be solved. When we study the competition system between product and product of competitor 1, we also apply the aforementioned approach to the case where phase separation of anhydrides is absent and the system is always homogeneous.

Obtaining reaction rate constants for the kinetic model
We study three succinate derivatives, which we label as P , C 1 and C 2 and their corresponding anhydrides. We fit the experimental measurements with the kinetic traces of each of the reaction cycles to determine the rate constants (Figs.3, 4 and 5).

Systems with one succinate derivative
Figure 3: First order deactivation and short lifetime of product of competitor 1. 50 mM competitor 1 fuelled with 100 mM EDC. The two curves, corresponding to the time trace of the product of competitor 1 and fuel concentration respectively, are globally fitted to obtain the reaction rate constants. The concentration profile of product of competitor 1 shows an exponential decay as it is not able to phase separate due to its high solubility of roughly 3000 mM. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. The three curves, corresponding to the time trace of the product, precursor and fuel concentration respectively, are globally fitted to obtain the reaction rate constants. The product concentration profile shows an exponential decay as it is not able to phase separate due to its high solubility of roughly 27 mM. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. The three curves, corresponding to the time trace of the product of competitor 2, competitor 2 and fuel concentration respectively, are globally fitted to obtain the reaction rate constants. The concentration profile of product of competitor 2 shows a linear decay as it is able to phase separate due to its low solubility of roughly 2 mM. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. The concentration of competitor 2 had to be re-adjusted to 57 mM in the theoretical kinetic model for the fitting procedure due to inaccuracies in the stock solution.
System with two competing succinate derivatives Figure 6: First order deactivation for the product and product of competitor 1. 50 mM precursor and 50 mM competitor 1 fuelled with 100 mM EDC. The three curves, corresponding to the time trace of two anhydrides and fuel concentration respectively, are globally fitted to obtain the reaction rate constants. The concentration profiles of product and product of competitor 1 both show exponential decay as neither is able to phase separate, and competition for fuel results in reduced yields for both anhydrides. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model.
Having obtained the reaction rate constants from the above fits in Fig. 6, we proceed to see how the traces change with changing competitor 1 (C 1 ) concentration and using fixed precursor (P ) and EDC concentrations (Fig.7).  respectively. 50 mM precursor (P ) and 50 mM competitor 2 (C 2 ) fuelled with 100 mM EDC. We globally fit five curves corresponding to the time trace of two anhydrides, two succinate derivatives and fuel concentration, respectively, to obtain the reaction rate constants. The two outlier points in the precursor concentration trace are omitted however. The linear decay of product of competitor 2 shows that it phase separates and the non-linear decay of the product implies that it partitions in the droplets, which we refer to as co-phase separation. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model.
Having obtained the reaction rate constants from the above fits in Fig. 8, we proceed to see how the traces change with changing competitor 2 (C 2 ) concentration and using fixed precursor (P ) and EDC concentrations (Fig.9).

Figure 9
: Competitor 2 reduces the yield but increases lifetime of the product. A) 5 mM C 2 , B) 10 mM C 2 , C) 25 mM C 2 , D) 75 mM C 2 , E) 100 mM C 2 , F) 125 mM C 2 and 50 mM P fuelled with 100 mM EDC. Increasing competitor 2 concentration reduces the maximum yield of the product, but prolongs its lifetime, allowing us to label competitor 2 as host and precursor as parasite. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. Figure 10: Reaction rate constants. Readout for the chemical reaction rate constants used in the kinetic model as obtained from the global fits.

Summary of reaction rate constants
We summarise the reaction rate constants obtained from fits as mentioned above (Figs.3,4,5,6 and 8). We observed that in general there are only little differences in the deactivation rate constants of the single precursor system and to the system where two components compete about fuel. However, a noticeable effect is the reduction in the activation pathway rate constant, primarily k f of competitor 1 (C 1 ) in the competition system with precursor (P ) and that of the precursor (P ) in the competition system with competitor 2 (C 2 ). For both single and competition studies we keep the solubilities of the product and product of competitor 2 unchanged to values, 2.01 mM and 27.8 mM, respectively. We thus assumed that competition for the fuel affects the availability of fuel for specific succinate derivatives. For periodic fueling studies, the lower value k d was used for the product. For the Fig. 5 we use the lower value of k d and for Fig.1F in the main text, we use the higher value of k d for product of competitor 2 as it was obtained from fitting routine.

Lifetime and half lifetime measurements for the product
The lifetime of the anhydrides is defined phenomenologically by threshold concentrations. The lifetime corresponds to the time period when the anhydride concentration is above this threshold. A threshold concentration can lose robustness against experimental measurements if it lies in the tailing regime of an exponential decay. It also becomes an injudicious choice if set to a high concentration value which is unachieved during the course of the experiment. Therefore, we choose threshold concentration values as described below such that the aforementioned difficulties are circumvented. The lifetime of product increases with increasing competitor 2 concentration. Markers represent lifetimes calculated using HPLC data; solid lines represent data calculated using the theoretical kinetic model.
We have performed two sets of competition experiments where precursor competes with competitor 1 and competitor 2, respectively. The important difference between the two competitors is that competitor 2 has a significantly lower solubility of about 2 mM compared to competitor 1 which has solubility of roughly 3000 mM, and therefore the former is capable of phase separation in the working concentration ranges. The lifetime measurements, corresponding to Fig. 3B in the main text, are done by tracking a threshold concentration of 2.01 mM for the product. Similarly, corresponding to Fig. 3E in the main text, the lifetime measurement of the product is done by tracking the threshold concentration 2.01 mM for the product of competitor 2, as it highlights the lifetime of the droplets. We believe that the lifetimes of both the host (product of competitor 2) and parasite (product) are the same in this case. However this form of measurement is undefined for small concentration values of competitor 2 since droplets are short-lived and the product survives after droplet dissolution. In these concentration ranges we track the threshold concentration for the product instead of the product of competitor 2. We therefore choose to measure the half lifetime of the product which has a robust definition in either competition system.
The half lifetime is measured as the time difference between the time points of the maximum value and the half maximum value of the product concentration, respectively. In the experimental data, we locate the maximum value of the product concentration, calculate its halved value and determine the time point by linearly interpolating between the immediate high and low value around the half maximum. Refer to Fig. 12 for visualisation. Half lifetime of product as a function of competitor 1 concentration. The half lifetime reduces due to competition for fuel in the system. B) Half lifetime of product as a function of competitor 2 concentration. The half lifetime increases due to protection in droplets from hydrolysis, but approaches saturation owing to a fixed fuel concentration of 100mM. The gray line denotes that in absence of co-phase separation in the competition system with competitor 2, the half lifetime would decrease solely due to competition for the fuel. Markers represent half lifetime calculated using HPLC data; solid lines represent data calculated using the theoretical kinetic model.

Composition of droplets
We quantified the percentage of the droplet material composed of the product for three different conditions with increasing competitor 2 concentration. Increasing concentration of the host reduces the maximum concentration of the parasite in the droplets. Also the maximum value of product monotonically decreases with time as the total volume of droplets keeps decreasing. Markers represent ratio calculated using HPLC data; solid lines represent data calculated using the theoretical kinetic model. Figure 14: Total droplet volume increases and total hydrolysis rate decreases in the dilute phase allowing longer survival of the parasite. A) The kinetic orbits corresponding to different concentrations of competitor 2. The kinetic orbit crosses fewer tie lines at maximal competitor 2 concentration (125 mM), suggesting the solubility change is not drastic and allows nearly zeroth order decay for both anhydrides. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. B) The total droplet volume increases with increasing the competitor 2 concentration and the product of competitor 2 starts acting like a host to protect the parasite inside the droplets from hydrolysis. C) The equilibrium concentration of the product in the aqueous phase, c out , decreases with increasing competitor 2 concentration. D) The hydrolysis rate of the product in the dilute phase decreases with increasing competitor 2 concentration. E) The hydrolysis rate of the product of competitor 2 in the aqueous phase sets the offset of the total hydrolysis rate, and it increases with increasing competitor 2 concentration. F) The total hydrolysis rate of both anhydrides outside the droplets. The stars in B, C, D, E and F denote the time-point of dissolution of droplets. Figure 15: Survival and consequent enrichment of both the product of competitor 2 (host) and the product (parasite) due to co-phase separation. A) The time traces of both anhydrides (host and parasite with initial concentrations 100 mM and 50 mM, respectively), for periodic fuelling at the rate of 60 mM every 30 minutes. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. Disagreement between the experimental data and the data calculated from the kinetic model for the host is more prominent after first cycle of fuelling due to formation of multiple droplets that settled at the bottom of the HPLC screw cap vial. B) The time traces of the host and the parasite under the experimental fuelling conditions to highlight the long time behaviour using the kinetic model. We observe the emergence of a non-equilibrium steady state for both the host and the parasite. C) The total volume of droplets also achieves a non-equilibrium steady state. D) The time trace representation in the phase diagram shows the oscillatory behaviour and approach to a non-equilibrium steady state. The inset depicts the limit cycle in the phase diagram achieved due to emergence of non-equilibrium steady states for both host and parasite (in the inset: last 3 cycles). Figure 16: The buffer capacity against fuel oscillations increases for the parasite with increasing host concentration. The buffer capacity defined as inverse of the deviation in concentration, i.e., (∆ −1 ) increases for the parasite with the increasing competitor 2 (host) concentration. The propensity to co-phase separate increases in the system with increasing the host concentration which allows for more protection of the parasite and thus less degradation. It leads to smaller deviations around mean psuedo-steady state concentration. The host's buffer capacity decreases, making it more susceptible to fluctuations. The blue and red solid line represent buffer capacity of host and parasite respectively. Initially the buffer capacity of the host is high due to lower mean concentration of the host and hence less deviation. As the mean concentration increases, the deviation around it also increases, thus reducing its buffer capacity. The opposite trend occurs for the parasite, upto 40 mM, following which the mean concentration of the product also increases due to the protection from the host droplets. The highlighted values of buffer capacity are corresponding to the data shown in the main text Figs. 4C and 4E, respectively.

Effects of activation rate constants on host-parasite identity
We have shown that the kinetic orbit of the average product concentrations in the phase diagram determines the lifetimes of the products and the composition inside droplets. The shape of the orbit is also affected by the rate constants. To illustrate this aspect we considered, precursor concentration of 50 mM and competitor 2 concentration of 100 mM fuelled with 100 mM EDC, as the experimental reference and swapped the rate constants related to the activation reaction pathway (Fig. 17).

Figure 17
: Host and parasite identity depends on the solubilities of the components and initial precursor concentrations. A) Kinetic orbits in the phase diagram corresponding to three different parameter sets in which we swapped the rate constants of the fuel-driven activation pathway (i.e., all rate constants except for the deactivation rate constants) and considered different concentrations of competitor 2 at a fixed EDC concentration of 100 mM and precursor concentration of 50 mM. Markers represent HPLC data; solid lines represent data calculated using the theoretical kinetic model. B) The ratio of product of precursor in droplets over time shows that the identity of host and parasite can change with swapped rates. Due to the higher solubility, the product is typically the parasite, except when it is at excess. In this case, the product starts as a host and transits to a parasite, as long as droplets do not dissolve beforehand. The star markers denote the time-point of droplet dissolution.