Biofunctionalized Materials Featuring Feedforward and Feedback Circuits Exemplified by the Detection of Botulinum Toxin A

Abstract Feedforward and feedback loops are key regulatory elements in cellular signaling and information processing. Synthetic biology exploits these elements for the design of molecular circuits that enable the reprogramming and control of specific cellular functions. These circuits serve as a basis for the engineering of complex cellular networks, opening the door for numerous medical and biotechnological applications. Here, a similar principle is applied. Feedforward and positive feedback circuits are incorporated into biohybrid polymer materials in order to develop signal‐sensing and signal‐processing devices. This concept is exemplified by the detection of the proteolytic activity of the botulinum neurotoxin A. To this aim, site‐specific proteases are incorporated into receiver, transmitter, and output materials, and their release, diffusion, and/or activation are wired according to a feedforward or a positive feedback circuit. The development of a quantitative mathematical model enables analysis and comparison of the performance of both systems. The flexible design could be easily adapted to detect other toxins or molecules of interest. Furthermore, cellular signaling or gene regulatory pathways could provide additional blueprints for the development of novel biohybrid circuits. Such information‐processing, material‐embedded biological circuits hold great promise for a variety of analytical, medical, or biotechnological applications.

. Design of the 3CPRO construct of module R.  Table S1. Parameters of the feedforward and the feedback system. 15 Table S2. Plasmids used in this study. 16 Table S3. Primers used for cloning of the plasmids in Table S2. 18 Suppl. Information S1. Mathematical Modeling

Derivation of the mathematical model
Mathematical models have been widely used to mechanistically understand and predict the performance of biochemical reaction networks. In the following we model the biomaterials based detection systems as biochemical reaction networks incorporating all causal dependencies of the network structure. We use non linear ordinary differential equations (ODEs) to describe the dynamic behavior of relevant system components termed states. Enzymatic reactions are modeled via Michaelis Menten kinetics, whereas basal activation and deactivation rates are implemented as mass action. The unknown model parameters are estimated from experimental data using a maximum likelihood approach. To determine parameter uncertainties in terms of parameter confidence intervals and thus the identifiability of the system we analyze the profile likelihood for each parameter. In the following, the two biomaterials based detection circuits describing a feedforward ( Figure 1A and B) and an integrated feedforward feedback ( Figure 1C and D) system are derived.
The model scheme of the feedforward system is shown in Figure S2A. For the sake of simplicity, anchors of released proteins are not shown. After the first contact with an active chain of botulinum toxin A, 3C protease is released from its polymer and triggers a linear reaction cascade finally resulting in the cleavage of one anchor of the output protein mCherry (v8, v10). Its second anchor is directly cleaved by the 3C protease (v9, v11) resulting in a feedforward loop which enhances the systems response. The model consists of 13 reactions with corresponding fluxes v1 v13 ( Figure S2B). Note that Michaelis Menten terms are re parameterized to facilitate fitting of enzymatic reactions in their linear range in relation to the substrate. The km x parameters we used for modeling therefore represent the inverse of the original Michaelis Menten constants Kmx. Given these fluxes the resulting set of ODEs can be formulated ( Figure  S2C) to describe the dynamic behavior of the system.
The structure of the integrated feedforward feedback detection circuit, referred to as feedback system, is outlined in Figure S3A. In contrast to the feedforward system, it is not the 3C protease but the active chain of botulinum toxin A that directly triggers the feedforward loop. Here, the signal is already integrated on the level of the TEV protease that is immobilized with two anchors, each containing cleavage sites for botulinum toxin A and Caspase 3. Moreover, a positive feedback loop is integrated in the system as the release of Caspase 3 is triggered by active TEV. This model consists of ten equations with corresponding fluxes v1 v10 ( Figure S3B), resulting in the ODE system depicted in Figure S3C. Note that the release of Caspase 3 and mCherry are both catalyzed by TEV protease and are modeled with the same rate constants for simplification.

Likelihood based parameter estimation and identifiability analysis
In order to build up a mathematical model accurately describing an observed biochemical network, unknown model parameters have to be determined. It is often difficult to directly measure these parameters, e.g., kinetic rates. However, they can be estimated from experimental data by maximum likelihood estimation.
The ODE system of our two models can be described with the following vectorized scheme . (2.1) Internal states are described as with initial concentrations and depend on unknown dynamic parameters , e.g. biochemical rates. Internal states are linked to experimental data by an observation function g , (2.2) which depends on observation parameters , e.g. offset and scaling parameters. In addition, we assume a constant Gaussian error ) (t   (0,  2 ) as measurement error, whose variance is estimated simultaneously with the dynamic parameters. The maximum likelihood estimation is performed as described in [1] in detail.
The resulting maximum likelihood estimate is further analyzed to assess uncertainty and identifiability of the optimal parameters. We therefor use the profile likelihood method [2] yielding confidence intervals for each parameter that are subsequently used for model reduction as described in [3] .

Implementation of single experiments
Four different experiments were used to calibrate the models. In the following their implementation is described and corresponding observation functions are defined. All initial states are set to zero if not stated otherwise.

Experiment 1: Characterization of the feedforward system
In this experiment the biomaterials based feedforward system outlined in Figure S2 was implemented to measure the relative amount of the output protein mCherry (OUT) under different concentrations of 3C protease, Caspase 3, TEV protease and botulinum toxin A. Thus, the initial concentrations of 3CPRO bound , Casp3 bound and TEV bound were set to the applied concentrations and initials of all other states to zero. The starting value of OUTbound was set to . (3.4) As observation function we used . (3.5) The measurement error was modeled with a constant Gaussian error with the standard deviation . Experimental data and model fits are shown in Figure S4A H.

Experiment 2: Characterization of the feedback system
In this experiment the biomaterials based feedback system outlined in Figure S3 was implemented to measure the relative amount of OUT again under different concentrations of 3C protease, Caspase 3, TEV protease and botulinum toxin A. Thus, the initial concentrations of 3CPRObound, Casp3bound and TEVbound were set to the applied concentrations and initials of all other states to zero. The starting value of OUTbound was set to . (3.5) As observation function we used . (3.6) The measurement error was modeled with a constant Gaussian error with the standard deviation . The experimental data and model fits are shown in Figure S5.

Experiment 3: Characterization of the transmitter module T2
In this experiment the T2 module of the biomaterials based feedback system was implemented to measure the release of TEV protease by different doses of Caspase 3 after four hours. The initial concentration of Casp3active was set to the different used concentrations. As applied in the experiment initTEV_bound was set to 500 RU and initial concentrations of 3CPRObound, Casp3bound and OUT bound were set to zero. The following observation function was used . (3.1) The measurement error was modeled with a constant Gaussian error with the standard deviation . Experimental data are shown in Figure 3E Experiment 4: Characterization of the output module O In this experiment the output module O of the biomaterials based feedback system was implemented to measure the release of OUT by different doses of TEV protease after 0, 6, 10 and 20 hours. The initial concentration of TEV was set to the different applied concentrations and initial concentrations of 3CPRO bound , Casp3 bound and Casp3 active were set to zero. The starting value of OUTbound was set to . (3. 2) The following observation function was used The measurement error was again modeled with a constant Gaussian error with the standard deviation . The experimental data are taken from Figure 3b in [1] .
Standard deviations as well as initial values for OUT bound differing from zero were estimated simultaneously with the dynamic parameters.

Fitting process and results
The feedforward and the feedback system have partly overlapping components and catalyzed reactions. Therefore, the two mathematical models share some of the parameters and fitting of the models to the experimental data was performed simultaneously. In total 29 parameters were fitted, including 21 dynamic parameters, two initial parameters, two scaling parameters and four error parameters. The set of estimated parameters is shown in Table S1.
Compilation, numerical integration, fitting and optimization of the model was performed with the MATLAB based free available software Data2Dynamics [4] . For the numerical integration of the ODE system we used the CVODES [5] solver and optimization was performed with the trust region based algorithm LSQNONLIN implemented in MATLAB [6] . The fitting was performed in logarithmic parameter space to scan the parameters over orders of magnitude. Moreover, multiple optimization runs with random sampled initial parameter sets were performed to ensure global convergence of the optimizer. Out of total 100 fits more than 80 % converged to the same lowest minimum indicating that the global optimum was found. The resulting model curves and data are shown in Figures S4 and S5. Shaded bands correspond to the estimated standard deviation of the Gaussian error model.
The identifiability analysis with the profile likelihood method showed that all inverse Michaelis Menten parameters (kmx), except for kmCasp3activation_3CPRO that was fixed to the value 1.70 10 +00 as determined in [1] , were practical non identifiability towards zero. As introduced in [3] this is a case for model reduction. Thus, affected km values were set to zero for the following analysis resulting in linear reactions. It turned out that all other parameters are identifiable ( Figure S6). 95 % point wise confidence intervals calculated with the profile likelihood method are listed in Table 1 reaching from  to  + .

Model validation
In order to test the model for robustness an eight fold cross validation was performed. The 64 different experimental conditions in the data, excluding the single module characterization (Exp. 1 and 2), were randomly distributed into eight groups of 8 conditions. Seven data groups were used for parameter estimation to predict the dynamics of the eighth group. All predicted trajectories were describing the unfitted data very accurately with almost no visible difference to the originally fitted curves. Estimated parameter values of these eight optimization runs with reduced number of data conditions are indicated by colored stars in Figure S6. More than 95 % of all estimates are within the 95 % confidence interval of the original estimate calculated with the full data set indicating a robust model structure.       Bacterial expression vector encoding the crosslinking OUT protein of the output hydrogel of the feedback system. [1,7] Crosslinker of madule O of the feedback system. pHJW4 PT7 His6 3CPRO Bacterial expression vector encoding His tagged human rhinovirus type 14 3C protease. [1,7] Characterization of 3CPRO mediated dissolution of module O of the feedforward system ( Figure 4C). unpublished Cloning of pHJW123.
pHJW123 PT7 His6 AviTag GyrB BTS EGFP 3CPRO The sequence of EGFP 3CPRO and one half of the backbone was amplified from pHJW183 using oligonucleotides oHJW5 and oHJW427. The second half of the backbone plus His6 AviTag GyrB BTS was amplified from pHJW81 using oligonucleotides oHJW6 and oHJW194. The two resulting fragments were assembled by Gibson cloning [8] , yielding a bacterial expression vector for the production of HRV14 3C protease fused to His tagged GyrB, SNAP25(141 206) and EGFP.
This work Supplementary Figure S1B. Cloning of pHJW178.
pHJW125 P T7 His 6 AviTag GyrB GyrB CCS BTS* TEV Q175N and Q178N mutations were introduced into BTS via site directed mutagenesis of pHJW14 using oligonucleotides oHJW486 and oHJW487. This resulted in a bacterial expression vector encoding TEV protease linked to tandem GyrB via a CCS and BTS* (SNAP25 (141 206, Q175N, Q178N)) containing linker.
This work Characterization of TEV mediated dissolution of material O of both systems ( Figure 4).

Plasmid Description Reference
Used for pHJW182 P T7 Casp3(3CPRO ind ) 3CS His 6 The TCS of pHJW181 was exchanged with 3CS by amplifying Casp3(3CPRO ind ) together with one half of the backbone using oligonucleotides oHJW6 and oHJW502, and the 3CS His tag containing linker region plus the second part of the backbone using oligonucleotides oHJW5 and oHJW501. The two resulting fragments were assembled by Gibson cloning, yielding a bacterial expression vector for Casp3 OFF of the feedforward system.