A Bayesian Approach to Extracting Kinetic Information from Artificial Enzymatic Networks

In order to create artificial enzymatic networks capable of increasingly complex behavior, an improved methodology in understanding and controlling the kinetics of these networks is needed. Here, we introduce a Bayesian analysis method allowing for the accurate inference of enzyme kinetic parameters and determination of most likely reaction mechanisms, by combining data from different experiments and network topologies in a single probabilistic analysis framework. This Bayesian approach explicitly allows us to continuously improve our parameter estimates and behavior predictions by iteratively adding new data to our models, while automatically taking into account uncertainties introduced by the experimental setups or the chemical processes in general. We demonstrate the potential of this approach by characterizing systems of enzymes compartmentalized in beads inside flow reactors. The methods we introduce here provide a new approach to the design of increasingly complex artificial enzymatic networks, making the design of such networks more efficient, and robust against the accumulation of experimental errors.


Materials and Instrumentation Materials
All chemicals and reagents were used as received from commercial suppliers without any further treatment unless stated otherwise.

Enzymes
The enzymes immobilized on PEBs were purchased from Merck Sigma-Aldrich. Specifically:

Flow setup
A custom made CSTR (see Figure S1 for the design schematic, effective volume 100 ) was charged with PEBs (in a ratio of 1:31 mg beads:injection volume). The openings of the reactor were sealed with Whatman Nuclepore Track-Etch polycarbonate membranes (5 poresize, cat. number 10417414) to prevent outflow of PEBs. To subject the CSTR to various flow conditions we used Cetoni Low-Pressure High-Precision Syringe Pumps neMESYS 290N and gastight Hamilton syringes, with all flowrates programmed using the Cetoni neMESYS software. For all methods of offline detection, the outflow of the CSTR was connected to a BioRad 2210 Fraction collector, collecting for 7.5 (absorbance) or 10 minutes (HPLC) per fraction. This experimental setup is shown in Figure S2a

Spectroscopy
Offline fluorescence Collected fractions were pipetted onto a microplate (Greiner Bio One, black, polystyrene, 96 flat bottom chimney wells) at 60 per well. Fluorescence measurements for enzyme activity determination were performed with a Tecan Spark 10M plate reader. The fluorescence intensity of wells containing 30-200 of the reaction mixture was monitored for 1-4 min (shaking 3s/orbital mode/amplitude 4mm) at 23°C using top or bottom reading mode, at / = 380nm/460nm for 7-amino-4-methylcoumarin (AMC) based substrates.
Offline absorbance Collected fractions were pipetted onto a microplate (Greiner Bio One, transparent, polystyrene, 96 flat bottom chimney wells) at 60 per well. The absorbance of each sample at 340 nm was measured in a TECAN SPARK M10 platereader. Using a calibration curve, the absorbance data was converted to NADH concentration.
Online absorbance Absorbance in flow experiments was continuously measured at the reactor's output with a custom made flow cell kindly provided to us by LabM8 (shown in Figure S2c), connected to an AvaLight 355 nm LED lamp. Absorbance at 380 nm was detected using an AvaSpec-2048 with 1,005 ms integration time and averaging at 200 times per recorded datapoint. Using a calibration curve, the absorbance data was converted to NADH concentration.

High-performance liquid chromatography
High-performance Liquid Chromatography (HPLC) was performed using Shimadzu NexeraX3/Prominence system under a 0.9 mL/min flow at 45°C with a Shimadzu WAX-1 column. For the Shimadzu system, an injection volume was subjected to a 25 min gradient program was used starting from 20 mM of potassium phosphate at pH 7.0 • 1 min -20 mM at pH 7.0 • 16 min -480 mM at pH 6.8 • 19 min -480 mM at pH 6.8 • 22 min -20 mM at pH 7.0 • 25 min -20mM at pH 7.0 Retention time of ADP is at 8.8 minutes, ATP at 15.5 minutes.

G6P-DH assay
Collected fractions were pipetted onto a microplate (Greiner Bio One, transparent, polystyrene, 96 flat bottom chimney wells) at 60 per well. To these wells was added an 20 of 5 mM of NADP+ and 1 of 500 u/mL of G6P-DH. The absorbance of each sample at 340 nm was continuously measured in a TECAN SPARK M10 platereader. Using a calibration curve, the absorbance data was converted to NADPH concentration. We assumed this NADPH concentration equivalent to G6P concentration.

Production and characterisation of polyacrylamide-enzyme beads
Microfluidic devices fabrication for beads production Figure S3: Schematic of the microfluidic device for droplets production (adaptated from 1 ) Microfluidic devices, designed according to the scheme shown in Figure S3 1 , were produced in silicon wafers by photo and soft lithography, with different orifice sizes (20-80 ) at the T-junction, for droplet formation. After that, the wafers were used for the production of polydimethylsiloxane (PDMS) devices. The PDMS replicas were separated from the wafer, and inlets and outlets of 1mm inner diameter were punched. After the replica was bonded by oxygen plasma treatment to a glass slide, channels were coated with 2% 1H,1H,2H,2H-Perfluoro-octlytriethoxysilane to make them hydrophobic 2 .

Method 1: Enzyme-first functionalization
This method was used for the synthesis of PEBs containing Trypsin.

Synthesis 6-acrylaminohexanoic acid succinate (AAH-Suc) linker
N-6-acryloyl amido hexanoic acid (30 mg, 0.172 mmol, 1 eq) and N-Hydroxy succinimide (20mg, 0.172 mmol, 1eq) were added in a large glass vial and brought under N2 atmosphere. Subsequently, they were dissolved in 1 mL of anhydrous DMF and cooled down to 0°C. Then, N, N'-dicyclohexylcarbodiimide (DCC, 38mg, 0.172mmol, 1eq) was added, and the vial was sealed and left to stir at 0°C for one hour. Subsequently, the reaction was brought to 4°C and stirred for 16h. The reaction mixture was assumed quantitative and immediately used for enzyme functionalisation without further analysis 3 .

Enzyme functionalization
The desired enzyme (1 eq) was added to a falcon tube containing 0.1 M NaHCO3 (4mL), followed by addition of AAH-Suc (7.5 eq) in DMF. The mixture was left to stir at 21°C for 1 h, after which the reaction mixture was dialysed and lyophilised.

General Enzyme Immobilization During Polymerization (EIDP) method
An emulsion of monodisperse water-in-oil droplets was produced by using the microfluidic device described above with a 20 -wide T-junction, following an adapted procedure from Rivello et al. [3] The reagents phase used for hydrogel formation was composed by 9.7% (w/v) acrylamide, 0.4% (w/v) bisacrylamide, 1.5% (w/v) 2,2'Azobis(2-methylpropionamidine) dihydrochloride (AAPH) and 0-8% (v/v) of a solution of functionalized enzyme (100 ). The oil phase was composed of fluorinated fluid HFE-7500 (3M) and 1.5% (v/v) Pico-Surf™ 1. The flow rates used to produce monodisperse beads of 50 average diameter were 600 /ℎ (Qw) for the reagents phase and 900 /ℎ for the oil phase (Qo). The emulsion of droplets created S6 in the microfluidic device was collected in an Eppendorf with 100 uL mineral oil to avoid evaporation and breaking the emulsion and polymerised for 10 minutes under UV light. After polymerisation, beads were washed three times with 20% (v/v) 1H,1H,2H,2H perfluorooctanoic (PFO) in HFE-7500 oil to break the emulsion. The obtained beads were then washed three times with 1% (v/v) Span 80 in heptane, three times with 0.1% (v/v) Triton X-100 in miliQ and three times with Milli-Q. After every washing step, the beads were centrifugated 30 seconds at 5000 rcf, and the supernatant was removed by pipetting. The resulting beads were freeze-dried and re-dissolved in miliQ at a concentration of 0.0322mg/ . Bead size was obtained after freeze-drying by imaging with a light microscope with a 40x objective lens. The average bead size was 50 .

Method 2:
This method was used for the synthesis of PEBs containing GDH, HK, and G6PDH.

Empty bead production method
The microfluidics device described above was used to produce gel beads with a 20 wide T-junction, following an adapted procedure from 2 . The gel solution phase consisted of 9.6% (w/v) acrylamide, 0.4% (w/v) N,N'-methylenebisacrylamide, 0.5% (w/v) acrylic acid and 1.5% (w/v) 2,2'Azobis(2-methylpropionamidine) dihydrochloride. The oil phase contained 1.5% (v/v) Pico-Surf™ 1 in fluorinated fluid HFE-7500 (3M). The flow rates for gel phase and oil phase were 600 /ℎ and 900 /ℎ, respectively. The outflow emulsion was collected in the tube which was filled with 100 mineral oil. Afterwards beads were polymerised using UV lamp for 10 minutes at 70% gain. After polymerisation the lowest layer containing fluorocarbon phase was carefully removed with a P200 pipette. The remaining beads were washed 3 times with 20% (v/v) 1H,1H,2H,2H-Perfluoro-1-octanol in HFE-7500 (3M), then 3 times with 1% (v/v) Span 80 in hexane, 3 times with 0.1% (v/v) Triton X-100 in Milli-Q and finally 3 times with Milli-Q. Every washing step was finalised with mixing the tube using vortex, centrifuging at 5000 x g for 3 min and removing the layer which was not containing beads. Furthermore, beads were flash frozen using nitrogen and freeze dried overnight. After re-wetting beads, their size was determined using light microscope with 40x magnitude objective. The average size was 50 in diameter.

General Enzyme Immobilisation after Polymerisation (EIAP) Method
Empty acrylamide beads were re-dissolved in Milli-Q at a concentration of 0.0322 mg/ . 1-(3-Dimethylaminopropyl)-3-ethylcarbodiimide hydrochloride (100 mM) and N-Hydroxysuccinimide (100 mM) were added to the reaction mixture. The total volume of the activation solution was 5-fold of the beads volume. The reaction mixture tube was put on the roller bank for 30 min. After this, beads were centrifuged at 5000 x g for 3 min. The supernatant was carefully removed using P200 pipette. Beads were washed 3 times by adding Milli-Q, mixing using vortex, centrifuging and removing the supernatant. Specific to each batch of PEBs, a certain amount of enzyme was added to the beads (see Table S1). The tube was put on the roller bank for 2 h coupling step. Sequentially, beads were washed 8 times by adding Milli-Q, centrifuging and removing the supernatant. Finally, beads were flash frozen using nitrogen and freeze dried overnight.

Overview of experiments
Processed datafiles from these experiments can be found on the accompanying github repository in the data folder as csv-files. The files are directly used during analysis performed in the Jupyter notebooks.
In the tables below, the volume refers to the volume of PEBs suspension injected into the CSTR, where the PEBs suspension itself is created by suspending 1 mg of dry PEBs in 31 buffer, or other quantities in the same ratio (1:31). The batch refers to the batch of PEBs, which can contain different effective enzyme concentation. See Table S1 for an overview of all PEB batches.

Overview of computational methods
Python scripts and Jupyter notebooks were used to create the Bayesian models and perform inference and predictive sampling. All computational studies were performed with Jupyter notebooks. Datasets were loaded in from csv-files with Pandas, and if relevant, concatenated together into larger objects.
The Bayesian model, including the determination of prior probabilities and likelihood function, was created using PyMC3 4 . Generally, prior probabilities for Michaelis-Menten parameters were chosen as uniform distributions over a specified interval. These distributions were used as uninformative priors to ensure no subjective information would enter the model, while garanteeing correct sampling and estimation of parameters. Priors for the uncertainty estimations (denoted by sigma in the notebooks) were given an exponential distribution, which also acted as an uninformative distribution, while garanteeing correct sampling and convergence. In larger models, where multiple likelihoods were combined, hyperpriors were placed on the and priors to increase convergence of the sampling algorithm.
All sampling was peformed using the No-U-Turn Sampler (NUTS) 5 , which is an adaptive step-size Hamiltonian Monte Carlo sampler. When (automatically, or via a custom operator) the gradients of the likelihoods with respect to the kinetic parameters were given, this sampling method is much more efficient then a classical Metropolis Monte Carlo sampler, showing faster convergence and requiring less samples for precise posterior estimations. Generally, sampling was performed using 4 or 8 independent chains on 4 or 8 cpu cores, all with 1000 tuning steps, and 1000 sampling steps, and a target step acceptance probability of 0.95. These values were found to yield good sampling results without becoming computatially inefficient.
For likelihood calculations of partially observable ERNs, a custom steady-state likelihood operator was written in Theano, according to description given in the manuscript. This operator is freely available from the public github repository at https://github.com/mgbaltussen/BayERN.
The samples obtained from the posterior distribution were further analysed using standard statistical tools in Python, the Numeric Python package and the Scientific Python package (NumPy and Scipy) 6 . To ensure the accessibility and reproducibility of these results the datasets and Jupyter notebooks, used for the analysis and creation of figures found in the publication, are made available as additional Supporting Information, and directly on github at https://github.com/huckgroup/Bayesian-enzymatic-networks_manuscript_2022. This repository also includes three example notebooks with additional explanation (example_simulation.ipynb, example_simple.ipynb, and example_complex.ipynb) and version information for all software dependencies.

Extended iterative combination of experiments
See Figure S4 for

Model details & sampling diagnostics
In the following section, the mathematical models on which the analysis is based are described in more detail, alongside additional information on sampling statistics and diagnostics. All three sections have a corresponding Jupyter notebook on the github repository at https://github.com/huckgroup/Bayesianenzymatic-networks_manuscript_2022.

Model
The cleavage of R-AMC to AMC by the peptidase Trypsin was modelled following a Michaelis-Menten rate equation with an additional term describing uncompetitive inhibition by the inhibitor AAA-AMC. This results in the following system of ODE's: As prior for the noise estimate sigma, an uninformative exponential distribution is used: The full posterior is described by where [ ] = ( , ) is the steady-state solution of f(C, , ). Table S4: Sampling statistics for manuscript Figure 1 Chains Tuning steps Draws Mean acceptance rate Divergences 0 4 1000 1000 0.948 0 S13  Figure S5: Diagnostic sampling trace plots of all kinetic parameters inferred in manuscript Figure 1 S15

Diagnostics
Combining diverse experimental datasets

Model
The system of GDH and HK PEBs was modeled as a competitive two-reaction system, where GDH catalyses the reaction G+NAD → GdL+NADH, and HK catalyses G+ATP → G6P+ADP. This system is described by the following set of ODE's:  The 0-value was not included in the domain of these priors because it is highly unlikely and causes problems with divergent sampling.
The turn-over number is not only different for both enzymes, but can also vary between batch of enzymes or PEB's. To somewhat constrain the inferred values and ensure that the sampler can properly converge, the priors were modelled as uniform distributions with an upper boundary sampled from a gamma-distribution.
As priors for the noise estimate sigmas for each experiment, an exponential distribution was used. To help the sampling procedure, an exponential hyperprior was used for the distribution parameters.   Comparing reaction mechanism hypotheses

Models
The system of G6PDH PEB's was modeled as the reaction G6p + NAD → G6PdL + NADH. Every hypothesis under investigation followed a different proposed rat equation, but the priors were kept as equivalent uniformative priors for every hypothesis: As prior for the noise estimate sigma, an uninformative exponential distribution was used:

Hypothesis comparison
Comparison of the predictive value of each hypothesis was done based on a PSIS-LOO cross-validiation calculation, with exact LOO cross-validation for observations where the Pareto shape exceeded 0.5. The exact LOO was performed by by direct resampling of the model with the selected observation left out.
Results of this procedure can be found in the table below.         Figure