A two-state photoconversion model predicts the spectral response dynamics of optogenetic systems

In optogenetics, light signals are used to control genetically engineered photoreceptors, and in turn manipulate biological pathways with unmatched precision. Recently, evolved photoreceptors with diverse in vitro-measured wavelength and intensity-dependent photoswitching properties have been repurposed for synthetic control of gene expression, proteolysis, and numerous other cellular processes. However, the relationship between the input light spectrum and in vivo photoreceptor response dynamics is poorly understood, restricting the utility of these optogenetic tools. Here, we advance a classic in vitro two-state photoreceptor model to reflect the in vivo environment, and combine it with simplified mathematical descriptions of signal transduction and output gene expression through our previously engineered green/red and red/far red photoreversible bacterial two-component systems (TCSs). Additionally, we leverage our recent open-source optical instrument to develop a workflow of spectral and dynamical characterization experiments to parameterize the model for both TCSs. To validate our approach, we challenge the model to predict experimental responses to a series of complex light signals very different from those used during parameterization. We find that the model generalizes remarkably well, predicting the results of all categories of experiments with high quantitative accuracy for both systems. Finally, we exploit this predictive power to program two simultaneous and independent dynamical gene expression signals in bacteria expressing both TCSs. This multiplexed gene expression programming approach will enable entirely new studies of how metabolic, signaling, and decision-making pathways integrate multiple gene expression signals. Additionally, our approach should be compatible with a wide range of optogenetic tools and model organisms. Significance statement Light-switchable signaling pathways (optogenetic tools) enable precision studies of how biochemical networks underlie cellular behaviors. We have developed a versatile mathematical model based on a two-state photoconversion mechanism that we have applied to the E. coli CcaSR and Cph8-OmpR optogenetic tools. This model enables accurate prediction of the gene expression response to virtually any light source or mixture of light sources. We express both optogenetic tools in the same cell and apply our model to program two simultaneous and independent gene expression signals in the same cell. This method can be used to study how biological pathways integrate multiple inputs and should be extensible to other optogenetic tools and host organisms.


TCS photoconversion model
We constructed a TCS 'sensing model' (Methods) by adding terms for production of new ground state photoreceptors (S g ) at rate k S and dilution of S g and active state photoreceptors (S a ) at rate k dil to the two-state model (Fig. 1a). The sensing model accepts any n light (λ) input and produces S g and S a populations as an output (Fig. 1b,c). The ratio S a /S g feeds into an 'output model' comprising a phenomenological description of TCS signaling and a standard model of output gene expression (Fig. 1c). The TCS signaling model (Methods) describes a pure time delay ( ) and Hill-function mapping ( (̂,̂, , )) between S a /S g and output gene production rate (k G ). In our initial experiments, we utilize superfolder GFP (G) as the output and quantify its expression level in Molecules of Equivalent Fluorescein (MEFL). ̂ is the range of possible k G values, ̂ is the minimum value of k G , n is the Hill parameter, and k is S a /S g ratio resulting in 50% maximal system response. Together, these terms capture SK autophosphorylation, phosphotransfer, RR dimerization, DNA binding, promoter activation, and G production. G is degraded in a first-order process with rate k dil (Methods), and has a minimum concentration =̂/ dil and concentration range =̂/ dil given a constant cell growth rate.

Light source model
Most light sources have a fixed spectral flux density (i.e. output spectrum) that scales with light intensity (I, µmol m -2 s -1 ). For such light sources, we can write light =̂l ight ⋅ where ̂l ight is the output spectrum at 1 µmol m -2 s -1 . To quantify the overlap between light and for a given photoreceptor, we introduce ̂ as the photoconversion rate per unit light intensity (min -1 [µmol m -2 s -1 ] -1 ). Then, for a given light source, = ⋅ ∫ ⋅̂l ight = ⋅̂. That is, k 1 and k 2 take on values proportional to light intensity.

Dynamical and spectral characterization of CcaSR
We designed a set of four gene expression characterization experiments (File S1-2, Note S1-2) to train the TCS photoconversion model for CcaSR (Fig. 2a). First, we quantify activation dynamics by preconditioning E. coli expressing CcaSR (Fig. S1) in the dark, introducing step increases in green light (centroid wavelength = 526 nm, Table S1-3, File S3-4, Note S3) to different intensities, and measuring sfGFP levels over time by flow cytometry (Methods, Fig.  2b, S2). Second, we measure de-activation dynamics by preconditioning the cells in different intensities of green light and measuring the response to step decreases to dark (Fig. 2c, S2). Third, we measure the ground state spectral response by exposing the bacteria to 23 LEDs with spanning 369 to 958 nm at over three orders of magnitude intensity (Methods, Fig. 2d, S3, Table S1-3, File S3-4, Note S3) and measuring sfGFP at steady state. Finally, we measure the activated state spectral response by repeating the previous experiment in the presence of a constant intensity of activating light (Fig. 2e, S3).

CcaSR model parameterization
We used nonlinear regression (Methods , Table S4) to fit the model to these data. While the resulting parameters recapitulate the known properties of the system ( Fig. 2f-g, S4), the value of the Hill parameter k is weakly determined ( Table S4). In particular, alterations in k from the best-fit value can be compensated for by changes in ̂1 and ̂2 (Fig. S5). Thus, we cannot confidently determine the absolute rates of forward and reverse photoconversion. Nonetheless, fixing k at its best fit value results in model predictions that quantitatively agree with the experimental measurements (Fig. 2b-e). However, the ultimate validation of this approach involves predicting the response of CcaSR to a wide range of spectral and dynamical light inputs different from those used in parameterization.

Spectral validation of the CcaSR photoconversion model
To predict the response of an optogenetic tool to a given light source, knowledge of is required. To estimate for CcaSR, we used non-linear regression to fit a cubic spline to the previously determined photoconversion rates for each of the 23 LEDs (Methods Fig. 3a, Fig.  S6-7). Importantly, our regression procedure considers the response of CcaSR to the full spectral output of each LED, not just its centroid wavelength. To validate the resulting estimate, we measured ̂l ight for a previously untested set of eight color-filtered white light LEDs designed to have complex spectral characteristics (Table S1-3, File S3-4, Note S3) and calculated an expected ̂ for each (Fig. 3b). In combination with the remaining model parameters (Fig. 2f), we used these ̂ to predict the steady state intensity dose-response to these eight LEDs in the presence and absence of activating light ( = 526 nm). These predictions are remarkably accurate for LEDs 1-5 (root-mean-square errors (RMSEs) from 0.11 to 0.18, Methods), which drive sfGFP to high levels, and 7 and 8, which drive low expression (RMSE = 0.14 and 0.18, respectively), but slightly less so for LED 6 (RMSE = 0.26), which drives sfGFP to an intermediate expression level (Fig. 3c).

Dynamic validation of the CcaSR photoconversion model
Our biological function generator method constitutes a rigorous validation of the predictive power of a model because the light inputs and gene expression outputs are temporally complex and cover a wide range of levels. To validate our CcaSR photoconversion model, we first designed a challenging reference gene expression signal (Fig. 4, File S5). The signal starts at b and then increases linearly (on a logarithmic scale) over 90% of the total CcaSR response range over 210 min. After a 60 min. hold, the signal decreases linearly to an intermediate expression level over another 210 min. Using this reference, we then used the model to computationally design four light time courses each with different LEDs or LED mixtures (Methods, File S6). "UV mono" utilizes a single UV LED ( = 389 nm) (Fig. 4a) to demonstrate control of CcaSR with an atypical light source. "Green mono" uses the = 526 nm LED (Fig. 4b) to demonstrate predictive control with a typical light source. "Red perturbation" combines "Green mono" with a strong red ( = 657 nm) sinusoidal signal (Fig. 4c) designed to demonstrate the perturbative effects of alternative light sources. Finally, in "Red compensation", the "Green mono" time course is re-optimized to compensate for the impact of "Red perturbation" (Fig. 4d, Methods).
The model accurately predicts the response of CcaSR to all four light signals (Fig. 4). "Mono UV" presents the greatest challenge, resulting in an RMSE of 0.15 (Fig. 4a). We suspect that prediction errors in this program are due to PCB photodegradation, as we observed no significant toxicity via bacterial growth rate, and the prediction remains accurate until UV reaches maximum intensity (20 µmol m -2 s -1 ). "Green mono" (Fig. 4b) results in the lowest error (RMSE = 0.038), which is expected because this LED was used to perform the dynamic calibrations (Fig. 2b,c). As intended, "Red perturbation" results in an enormous deviation from the reference signal (Fig. 4c), and the model accurately predicts this effect (RMSE = 0.081). Finally, "Red compensation" demonstrates that the effect of the perturbation can be eliminated using our model (Fig. 4d, RMSE = 0.078).

Cph8-OmpR photoconversion model
To evaluate the generality of our approach, we repeated the entire workflow for Cph8-OmpR ( Fig. S8-13, Table S5, File S6-7). Though CcaSR and Cph8-OmpR are both photoreversible TCSs, they have different photosensory domains, ground state activities, and dynamics. To account for the fact that Cph8-OmpR is produced in an active ground state, we used a repressing Hill function (Methods). The model again fits exceptionally well to the experimental data (Fig S8-11). Unlike CcaSR, which exhibited no detectable dark reversion (Fig. 2f), Cph8-OmpR appears to revert in 1/2 = ln 2 / dr = 5.5 min (Fig. S8f). As before, k is underdetermined, and we chose the best-fit value ( Table S5). The Cph8-OmpR model performs similarly to its CcaSR counterpart in the spectral validation experiments (Fig. S12), and demonstrates greater predictive control in the dynamical validation experiments (Fig. S13).

Development of a CcaSR, Cph8-OmpR dual-system model
We engineered a three-plasmid system (Fig. S1) to express CcaSR and Cph8-OmpR in the same cell with sfGFP and mCherry outputs, respectively (Fig. 5a). Because the photoconversion parameters are a property of the photoreceptors themselves, we left them unchanged. To recalibrate for mCherry (quantified in Molecules of Equivalent Cy5 (MECY)) and any changes due to the new cellular context, we measured the steady state levels of the sfGFP and mCherry at different combinations of green ( = 526) and red ( = 657) light ( Fig.  5b, S14, File S8) and refit the Hill function parameters ( Table S6). The dual-system model accurately captures the experimental observations from the characterization dataset (Fig. 5b).
To validate the dual-system model, we again used the biological function generator approach (Fig. 6). We designed a series of four dual sfGFP/mCherry expression programs to increasingly challenge the model: "Green mono" using only green light and intended only to control CcaSR (Fig. 6a), "Red mono" using only red light and intended to control only Cph8-OmpR (Fig. 6b), "Sum", a simple combination of the first two programs (Fig. 6c), and "Compensated sum" where the green light time course is re-optimized to account for the presence of the red signal (Fig. 6d) as before (Methods). Due to the minimal response of dualsystem Cph8-OmpR to green light (Fig. 5b), there was no need to adjust the red program to compensate for the presence of green light. The validation experimental results (Fig. 6) show that our dual-system model accurately captures both the sfGFP and mCherry expression dynamics. The CcaSR predictions are nearly as accurate as the single-system experiments ( Fig.  4), and the Cph8-OmpR results match single-system accuracy ( Fig. S12-13), demonstrating the extensibility of our approach to multiple optogenetic tools.

Multiplexed biological function generation
Finally, we designed and experimentally implemented four multiplexed sfGFP/mCherry expression functions representing classes of signals useful for gene circuit characterization (File S5-6). "Dual sines" illustrates that two gene expression sinusoids with different offsets, amplitudes, and periods can be composed without interference (Fig. 7a). Variations of this combination of signals could be used to perform frequency analysis of multiple nodes in a gene network. "Sine and stairs" demonstrates that our approach can generate two completely different gene expression signals at the same time (Fig. 7b). "Dual stairs" demonstrates that the ratio of two proteins can be varied over a remarkably wide range (Fig. 7c). Finally, "Time-shifted waveform" (Fig. 7d) demonstrates that our approach can be used to characterize genetic circuits where time-delays are critical, such as those involved in cellular decision-making.

Discussion
In this study, we demonstrate the first use of a mechanistic model of wavelengthdependent photoconversion to characterize and control light responsive signaling pathways in vivo. Additionally, we develop a standard set of characterization and validation experiments to parameterize the model and demonstrate that it accurately predicts the spectral and dynamical performance of these optogenetic tools. We demonstrate that the models can be used with virtually any light source or mixture of light sources as long as their emission spectra are known. Finally, we exploit this unique predictive power to demonstrate the first programming of two independent gene expression signals by accounting for inherent cross-talk in the action spectra of the two optogenetic tools that would otherwise impede such efforts.
Our TCS photoconversion model is superior to current alternatives by several key criteria. First, like our previous model 11 , it is quantitatively predictive and requires no parameter recalibrations from day-to-day. However, while our previous model is restricted to a single light source, our current model generalizes to virtually any light source or mixture of light sources. Second, our TCS photoconversion model is compatible with photoreceptors with very different action spectra, opposite ground versus active state signaling logic, and dramatically different dark reversion timescales. Third, our current model modularly decouples the processes of sensing (photoconversion) and output (signal transduction and gene expression). The sensing model component (Fig. 1a) should be compatible with a wide range of photoreceptors, including those in other organisms, because the core two-state photoswitching mechanism is used to describe their performance in vitro. Then, to describe optogenetic tools based upon those photoreceptors, our TCS output model can be replaced with alternatives appropriate to other pathways, as needed.
A major current problem in optogenetics is that tools developed in different studies are characterized using different culturing conditions, experiments, light sources, reporters, metrics, and so on. This lack of standardization makes it challenging to compare the performance features of different optogenetic tools on even a qualitative basis. The modeling and characterization approach we develop here could be used to make data sheets that describe the behavior of diverse optogenetic tools in standard units. This would enable researchers to choose the most appropriate tool for different applications. Additionally, shortcomings of specific tools could be identified, informing efforts to optimize performance by rational approaches such as protein design [16][17][18] .
Our approach should enable better control of optogenetic tools with alternative or highly constrained optical hardware used in many research laboratories. For example, many groups perform single cell optogenetic studies using fluorescence microscopes with severely restricted optical configurations. Alternatively, consumer projectors or tablet displays are potentially powerful, low cost hardware options for optogenetics 19,20 . The output spectrum of the light source can be measured and integrated into our workflow. After a simple recalibration (e.g. Fig.  5) to account for any changes due to the new growth environment, one should be able to predict and control the optogenetic tool using the new light source.
Oftentimes, it is desirable to simultaneously control an optogenetic tool while imaging a cell of interest using white light sources and excitation light for fluorescent reporters. Such alternative sources of illumination can have deleterious effects on the ability to control the optogenetic tool. However, if the nature of the alternative light signal is known, our approach can compensate for such perturbations (e.g. Fig. 6, 7). In silico feedback control has also been used to drive desired gene expression dynamics in optogenetic experiments [21][22][23] . The major benefit of this approach is that perturbations of unknown origin can be compensated by monitoring deviations in the output of an optogenetic tool relative to a reference. Our model is compatible with in silico feedback control.
While basic multichromatic control of optogenetic tools has been previously demonstrated 8,24 , the multiplexed biological function generation approach demonstrated here dramatically extends the capabilities of these systems, enabling implementation of several classes of experiments. First, the two-dimensional response of a genetic circuit or signaling pathway could be rapidly evaluated with high reproducibility and precision. For example, one could map the response of 2-input transcriptional logic gates 25 , which integrate the expression levels of two different transcription factors by systematically and independently varying their expression levels while measuring the gate output with a reporter gene. The dynamics of such gates are otherwise difficult to evaluate and seldom characterized 26 . Second, the input/output dynamics of a transcriptional circuit could be characterized as a function of the state of the circuit itself. For example, one could evaluate how well a synthetic transcriptional oscillator can be entrained 27,28 as a function of the strength of a feedback node. In this case, one optogenetic tool could be used for the entrainment, while the second was used to alter expression level of a circuit transcription factor regulating feedback strength. Third, transcription and proteolysis 29 could be independently controlled with two different optogenetic tools to alternatively program rapid increases or decreases in expression level. Such an approach could accelerate the gene expression signals that we have generated in this and our previous study 11 , enabling characterization of gene circuit dynamics on faster timescales. Finally, multiplexed biological function generation could be used to evaluate how the timing of expression of two genes impacts cellular decision making [30][31][32] . For example, in B. subtilis, the gene circuits that regulate sporulation and competence compete via a 'molecular race' in the levels of the corresponding master regulators 30 . By placing them under independent optogenetic control, the means by which their dynamics impact these cellular decisions could be evaluated more easily and rigorously.
Bacterial growth and light exposure Cell culturing and harvesting protocols were developed to ensure a high degree of precision and reproducibility in experiments both from well-to-well and from day-to-day (Note S1). Cells were grown at 37°C and shaken at 250 rpm throughout the experiment (Sheldon Manufacturing Inc. SI9R) with temperature calibrated and logged by placing a thermometer probe in a sealed 125 mL water-filled flask (Traceable Excursion-Trac 6433). Cultures were grown in M9 media supplemented with 0.2% casamino acids, 0.4% glucose, and appropriate antibiotics. Precultures were prepared in advance by freezing 100 µL aliquots of early exponential phase cultures (OD 600 = 0.1-0.2) grown in the same media conditions at -80°C (Note S2). Cultures were inoculated at low densities (typically OD600 = 1 × 10 −5 ) to ensure that final densities did not reach stationary phase (OD600 < 0.2). For each experiment, 192 cultures were grown in 500 µL volumes within 24-well plates (ArcticWhite AWLS-303008), sealed with adhesive foil (VWR 60941-126).
Experiments were performed using eight 24-well Light Plate Apparatus (LPA) instruments 34 , enabling precise control of two LEDs to define the optical environment of 192 cultures at a time. LPA program files were generated using Iris 34 and Python scripts.
LED measurement All LEDs were measured and calibrated (Note S3) using a spectrometer (StellarNet UVN-SR-25 LT16) with NIST-traceable factory calibrations performed on both its wavelength and intensity axes immediately prior to use for this study. A six-inch integrating sphere (StellarNet IS6) was used, enabling measurement of the total power output of each LED (in µmol s -1 ). The spectrophotometer was blanked by a measurement of a dark sample before each LED measurement. Measurements were saved as .IRR files, which contain the complete LED spectral power density light ( ) (µmol s -1 nm -1 ) in 0.5 nm increments as well as all setup parameters for the measurement (i.e. integration time and number of scans to average). These files were processed by Python scripts to calculate the LED characteristics, including the peak, centroid, FWHM, and total power. For spectral validation experiments, cinematic lighting filters (Roscolux) were cut, formed into LED-shaped caps, and fitted atop white LEDs ( Table S1).

Calculation of n light
Because the LEDs we utilize have fixed spectral characteristics, the spectral flux density (µmol m -2 s -1 nm -1 ) incident on the photoreceptors can be parameterized by the LED intensity (µmol m -2 s -1 ). The cultures are shaken throughout the experiment, and we assume that the cells are well mixed within the culture volume. Thus, the mean light intensity within the culture volume, light ( ), can be calculated by integrating the intensity throughout the volume of the well. Under the assumption of negligible light absorption by the culture sample (the M9 media is transparent, and the cultures are harvested at low density), this integral simplifies to become the total power of the LED (µmol s -1 ) divided by the cross-sectional area of the well. Given a well radius of 7.5 mm, we calculate light ( ) = light ( ) (7.5×10 −3 m) 2 ≈ 5.659 × 10 3 m −2 × light ( ).

LED calibration
Each of the approximately 700 individual LEDs used in the study were measured (Note S3), enabling compensation for variation in LED and LPA manufacturing (Table S1-3). Each LED was calibrated while powered from the same LPA socket used in experiments. First, a sample of LEDs were measured to identify the electrical current required to achieve an appropriate level of total flux, ∫ light ( )dλ. The amount of current required varied depending on the wavelength and manufacturer. The current was adjusted using the LPA 'dot-correction (DC)' to achieve a total flux approximately 20% above 20 µmol m -2 s -1 when the LED was fully illuminated. The appropriate DC level was determined for each LED model. Using these DC levels, the complete set of LEDs were measured. LEDs that produced a total flux below 20 µmol m -2 s -1 were remeasured at a higher DC level. This set of LED measurements was used to convert the desired intensity time course of each LED into a series of 12-bit grayscale values (i.e. 0-4095) used by the LPA. The LPA reads the grayscale values to produce the appropriate pulse-width-modulated (PWM) signal to achieve the desired intensities.
Bacterial sample harvesting Cultures were harvested for measurement (Note S1) after precisely 8 h growth by placing the 24well plates into ice-water baths. Each culture was then subjected to both an absorbance measurement to ensure consistent well-to-well and day-to-day growth, and flow cytometry for quantification of sfGFP or mCherry expression. Absorbance measurements were performed in black-walled, clear-bottomed 96-well plates (VWR 82050-748) in a plate reader (Tecan Infinite M200 Pro). Before fluorescence measurements were performed, culture samples were processed via a fluorescence maturation protocol to ensure measurements were representative of the total amount of produced fluorescent reporter 11 . Rifampicin (Tokyo Chemical Industry R0079), was dissolved in Phosphate-buffered saline (PBS, VWR 72060-035) at 500 µg/mL and used to inhibit sfGFP production during maturation.
Flow cytometry Population distributions of fluorescence were measured for each culture on a flow cytometer as previously described 11 . A calibration bead sample (Spherotech RCP-30-5A) in PBS was measured immediately prior to the culture samples from each experimental trial. At least 5,000 events were collected for the calibration bead sample, and at least 20,000 events were collected for each culture sample.
Flow cytometry data analysis Single-cell distributions of sfGFP fluorescence were gated, analyzed, and calibrated into MEFL and MECY units using FlowCal 35 . Measurements were gated on the FSC and SSC channels using a gate fraction of 0.3 for calibration beads and 0.8 for cellular samples 35 . Reported culture fluorescence values are the arithmetic means of the cellular populations.

Sensing model
The light sensing model can be described by the following system of ODEs: where the variables and rates have been described in the text and figures. Note that 1 and 2 are implicitly dependent upon time, as they are functions of the time-varying light environment of the sensors.
If we substitute for the fraction of active sensors, ≡ + , the system reduces to: where tot ≡ 1 + 2 + dil + dr .
This ODE can be solved analytically for a step-change in light from one environment to another. If the step-change occurs at time = 0, then 1 , 2 , and tot are all fixed for > 0. Given an initial sensor fraction (0) = 0 , we find.
This solution represents an exponential transition from an initial sensor fraction of 0 to a final fraction given by 1 tot with a time constant set by tot . As a result, we anticipate that the transition dynamics of ( ) will be slowest under zero illumination when tot = dil + dr . We also expect that the transition rates will be unbounded as intensity increases.
Finally, for multiple light sources, we simply linearly combine the photoconversion rates from each source: = , 1 + , 2 .
TCS signaling model We utilize a highly simplified model of TCS signaling and gene regulation. This model relates the production rate of the output gene ( ) to the active ratio of light sensors Generation of model simulations Simulations were produced by numerically integrating the system of ODEs using Python's scipy.integrate.ode method using the 'zvode' integrator with a maximum of 3000 steps.

Model parameterization
The CcaSR and Cph8-OmpR models were parameterized using global fits of the model parameters to the complete training data sets (Fig. 2b-e, Fig. S8b-e). The 'lmfit' Python package, which is based on the Levenberg-Marquardt minimization algorithm, was used to perform the fits and analyze the resulting parameter sets 36 . The fits were performed by minimizing the sum of the square of the relative error between each measured data point and the same point in a corresponding model simulation. Thus the form of the error metric utilized was error = ∑ ( data − model data ) 2 across the complete set of data points { data }.
Estimation of PCSs PCS estimates est. ( ) were constructed by linearly regressing a cubic spline to the experimentally determined photoconversion rates in order to produce a continuous PCS (Fig.  S6). The est. ( ) were produced by minimizing the error between unit experimental photoconversion rates ̂e xpt. (Fig. 2f, Fig. S8f) and spline-derived estimates ̂e st. = ∫ est. ⋅ light . The splines were constructed by establishing a series of integral constraints for the photoconversion rates, continuity constraints for the spline knots, and boundary constraints. As this problem contains more constraints than parameters, optimization is required. We used weighted least-squares with Lagrange multipliers to optimize each spline. To avoid overparameterization of the est. ( ), we used "Leave-one-out cross-validation (LOOCV)" to evaluate the performance of optimal splines with between 5 and 20 knots in order to determine the ideal number required for each est. ( ) (Fig. S7). Light program generator (LPG) algorithm The light program generator was used as previously described 11 . The only modification was to use simulations generated by the model described herein rather than the previous model.  The two-state photoreceptor model, which includes ground (S g ) and active (S a ) state photoreceptors (aka sensors), photoconversion rates k 1 and k 2 , and dark reversion rate k dr , is converted to a "sensing model" for in vivo environments by adding a S g production rate k S that captures both gene expression and holo-protein formation, and a dilution rate k dil due to cell growth and sensor degradation (Methods). The hollow blue pentagon represents a chromophore in the ground state while the filled blue pentagon represents that in the activated state. (b) Photoconversion rates are determined by the overlap integral of the spectral flux density of the light source (n light ) and the S g and S a photoconversion cross sections σ g and σ a (Methods). (c) The sensing model converts n light into the active ratio of light sensors S a /S g which feeds into an "output model" with a simplified model of TCS signaling that regulates the production rate k G of the target protein G, which is diluted due to cell growth and proteolysis as k dil (Methods). Training data for the full CcaSR system model (Fig. 1c). Experimental observations ("Expt.") and simulations of the best-fit model ("Model") are shown for each set. In particular, the response dynamics to step (b) increases from dark to eight different intensities and (c) decreases from eight different intensities to dark were evaluated using the = 526 nm LED. Time points are distributed unevenly to increase resolution of the initial response. (c-d) Steady-state intensity dose-response to a set of 23 "spectral LEDs" with spanning 369 nm to 958 nm. (c) Forward photoconversion is primarily determined by the response to the spectral LEDs. (d) Reverse photoconversion is analyzed by including light from a second, activating LED ( = 526 nm at 1.25 µmol m -2 s -1 ). The = 369 nm LED is not capable of reaching the brightest intensities, and thus those data points are not included. Light intensities are shown in units of 0.1 x log 2 µmol m -2 s -1 scale (e.g. a value of 1 corresponds to 10 x 2 1 = 20 µmol m -2 s -1 ). sfGFP fluorescence is calibrated to MEFL units (Methods). Each row of measurements in panels b-e was collected in a single 24-well plate. The 40 plates required to produce the training dataset were randomly distributed across eight LPAs over five separate trials (Methods, File S1-2). Each color patch represents the arithmetic mean of a single population of cells. (f-g) Best-fit model parameters produced via nonlinear regression of the model to training data (Methods , Table S4). ̂ are unit photoconversion rates (i.e. = ⋅̂, where I is the LED intensity in µmol m -2 s -1 ). Uncertainty in the least-significant digits are indicated in parenthesis.

Estimation of the CcaS photoconversion cross section and spectral validation of the CcaSR model. (a)
We estimate the continuous ground and active state PCSs of CcaS ( est. , lines) by regressing cubic splines to minimize the difference between experimentally-determined photoconversion rates (points) and those calculated via ̂e st. = ∫ est. ⋅ light (Methods, Fig.  S6,7). Error bars correspond to the standard error of the fits for the experimentally-determined photoconversion rates. The normalized spectral flux densities of the spectral LEDs are shown at bottom. (b) Using est. to predict photoconversion rates for light sources not in the spectral LED training set. Predicted photoconversion rates are integrated into the CcaSR model by keeping all other parameters (Fig. 2f) fixed, enabling prediction of the intensity dose-response of CcaSR to the new light source (i.e. ( ) . ). (c) Spectral validation of the CcaSR model and est. consists of prediction of the intensity dose-response for eight challenging, broad-spectrum light sources constructed by applying colored filters over white-light LEDs (Methods, Table S1-3, File S1-4). Measured n light , predicted ̂, measured and predicted intensity dose-response curves, and RMSE between model and prediction are shown for each LED (Methods). The forward and reverse intensity responses are determined using the filtered LED alone (circles) and in the presence of a second activating LED ( = 526 nm at 1.25 µmol m -2 s -1 , triangles). The simulated responses are determined using the calculated photoconversion rates (Methods). RMSE relative errors are expressed in log 10 decades (Methods). Data was collected across four LPAs, and the forward (circles) and reverse (triangles) intensity responses were collected over two separate experimental trials (Methods, File S1-2). Each data point represents the arithmetic mean of a single population of cells.  (File S5-6). (a) "UV mono". The LPG-generated UV light signal drives the CcaSR system along a trajectory predicted to follow the reference signal. (b) "Green mono". The green LED alone provides an optimized input signal. (c) "Red perturbation". The green LED provides the "Green mono" signal, while the red LED generates a sinusoidal perturbative signal (center) with a 240-minute period and 20 µmol m -2 s -1 peak-to-peak amplitude. (d) "Red compensation". The red perturbative signal is again present. However, the LPG redesigns the green light signal to account for its presence. Light signals are shown in units of log 10 µmol m -2 s -1 , and RMSE relative errors are expressed in log 10 decades (Methods). Error bars correspond to the standard deviation in fluorescence measurements over three independent experimental trials ( Table S4).  Fig. 2a (b) Training data for the multiplexed model ("Experiment", File S8) consists of a two-dimensional steady-state intensity dose-response to green ( = 526 nm) and red ( = 657 nm) light. The light intensities are logarithmically distributed, with the green light varying on a 0.05 x log 2 µmol m -2 s -1 scale (e.g. a value of -1 corresponds to 20 x 2 -1 = 10 µmol m -2 s -1 ) and the red light varying over a 0.05 x log 3 µmol m -2 s -1 scale (e.g. a value of -1 corresponds to 20 x 3 -1 = 6.67 µmol m -2 s -1 ). The different intensity ranges are used to maintain a high-resolution measurement despite the differences in the intensity dose-responses of the two systems. The four missing intensity values (white boxes) were not collected. The training data was used to re-fit the a, b, n, and k Hill function parameters for the CcaSR and Cph8-OmpR models (Table S6). Simulated steady-state responses to the same light environments for the best-fit dual-system models ( Table  S6) are shown ("Model"). mCherry fluorescence is calibrated to MECY units (Molecules of Equivalent Cy5, Methods). RMSE relative errors are expressed in log 10 decades (Methods). Data was collected in one experimental trial, and the 192 samples were randomly distributed across eight LPAs (Methods , Table S6). Each color patch represents the arithmetic mean of a single population of cells.  (Fig.5a) to time-varying signals of green ( = 526 nm) and red ( = 657 nm) light are compared to experimental results. Reference signals, light programs, and experimental data are as in Fig 4. (a) "Green mono". The green LED alone provides an optimized input signal for CcaSR. (b) "Red mono". The red LED alone provides an optimized input for Cph8-OmpR. (c) "Sum". The "Green mono" and "Red mono" programs are used simultaneously without any compensation, leading to a substantial deviation of the CcaSR output from the reference trajectory. (d) "Compensated sum". The "Red mono" program is used; however, the green light program is produced while incorporating red light program into the LPG (above). RMSE relative errors are expressed in log 10 decades (Methods). Error bars correspond to the standard deviation in fluorescence measurements over three separate experimental trials (Table S6).

Fig 7 -
Multiplexed biological function generation. The LPG is used to program CcaSR and Cph8-OmpR outputs to independently follow different reference signals. Red light ( = 657 nm) programs are optimized first using the LPG, and then the "Compensated" approach ( Fig. 6d) is utilized to generate the green light ( = 526 nm) program (Methods). The LPG, reference signals, light programs, and results are as in Fig 5d. (a) "Dual-sines". The sfGFP and mCherry reference trajectories are sinusoids with different periods, amplitudes, and offsets. (b) "Sine and stairs". The mCherry signal follows the same sinusoid in "Dual-sines," but the sfGFP reference is a stepped trajectory with several plateaus and increasing linear ramps. (c) "Dual-stairs". The sfGFP signal follows the same stair-shape in "Sine and stairs," however the mCherry response is a decreasing stair-shape. (d) "Time-shifted waveform". The sfGFP and mCherry reference trajectories both follow the same arbitrary waveform consisting of ramps, holds, and a sinusoid, with sfGFP trailing mCherry by 40-minutes. RMSE relative errors are expressed in log 10 decades (Methods). Error bars correspond to the standard deviation in fluorescence measurements over three independent experimental trials ( Table S6).

SI Materials and Methods
Detailed bacterial growth and light exposure protocol 1. Experiment initialization. This an approximately 2-hour process.
(a) Prepare program files on SD cards and load the SD cards into LPAs.
(b) Prepare 24-well plates for use (ArcticWhite AWLS-303008): Soak plates for at least 15 minutes in 70% EtOH. Triple-rinse the plates using DI water. Place rinsed plates onto a sheet of foil which will later wrap around the plate. Dry plates near a lit burner until the interior of the wells are dry (approximately 45 min). Use a paper wipe to dry the bottom of the plate and the foil. Wrap the plate in the foil.
(d) Distribute inoculated media into 24-well plates: Prepare a row of 1mL tips in a box so that only every-other tip is present (i.e. 6 tips total). Shake to homogenize inoculated media. Arrange the 8x 24-well plates with the wells open and uncovered. Pour more than half of the media into a 50mL disposable multichannel tray. Use a 1mL 12-channel pipettor with the previouslyprepared row of 6 tips to transfer 500 µL of media into each well. Pour remaining media into the tray and continue. Cover each plate with an adhesive-backed foil seal (VWR 60941-126).
(e) Load the LPAs and start the programs: Carry the 8x plates, 8x LPA lids, and 32x LPA wing nuts to the incubator containing the LPAs. (Note: a small carboard box or plastic container is helpful). Load the plates onto each LPA, ensuring that the plate is oriented correctly and is completely engaged with the LPA plate adapter. Place the lid onto each LPA, ensuring that the lid is oriented correctly. Engage 4x LPA wing nuts onto each device. Tighten the nuts evenly until the pressure of the rubber gaskets being compressed is felt (approximately 1-2 full turns after the nut engages with the lid). Start an 8-hour timer while synchronously releasing the reset button on one of the LPAs. Every 10 seconds, release the reset button on another LPA. Make sure to recall the order in which the LPAs were reset. This staggering of the start times enables the plates to be removed immediately when each of their programs ends. Allow the cells to grow in the LPAs at 37 °C for 8 hours.
2. Experiment completion and data collection. This is an approximately 4-hour process.
(a) 30 minutes prior to the end of the LPA programs, prepare a PBS+rifampicin solution: Fill a 250 mL beaker with 125 mL of a PBS solution in the pH7-7.2 range (VWR 72060-035). Add a stir-bar to the beaker and place on a stir-plate. Weigh 62.5 mg of rifampicin (rif, Tokyo Chemical Industry, R0079). Adjust the beaker to be centered on the stir plate. Adjust the rate of the stir plate to produce a vortex which doesn't quite reach down to the bar. Add the rifampicin to the middle of the vortex. Slide the beaker to be slightly off-center on the stir plate. This often encourages the vortex to lower, leading to the stir-bar clipping the vortex with each rotation. This Table S3. LED calibration results. Statistics generated during the calibration of the LEDs used in this study are shown. These results were calculated after the first-pass of calibration and thus are representative of the variation present in each set of LEDs given the same settings in the LPA (i.e. all LEDs used in each row of the table had the same "dot-correction" and "gray-scale calibration" value).

SI Files
File S1 Plate-by-plate overview of CcaSR experiments. This Excel spreadsheet contains a description of the light programs used for each 24-well plate used for the CcaSR experiments. On the "experiment" sheet, the "Run ID" and "Run name" columns contain unique identifiers for each experimental trial. The "Strain ID" is an identifier for the CcaSR strain. "LPA" is the name of the LPA device used. "Top LEDs" and "Bottom LEDs" correspond to the LED set identifiers used in the LED calibration archives. "Randomize" indicates whether the light program time points were randomized throughout the wells of the plate. "Program type" is an identifier corresponding to the type of light program ("aas" is a forward action spectrum, "ras" is a reverse action spectrum, "dta" is an activating step change, "atd" is a deactivating step change, "dv" is a dynamic validation, and "sv" is a spectral validation). The parameter columns further specify the details of the light program (intensities are in µmol m -2 s -1 , "Model" refers to an identifier for the system model used to generate the light program, reference signals and perturbation signals are identifiers corresponding to reference programs listed in File S3). For dynamic validation experiments, the model, reference signal, and perturbation signal (if present) are used to generate light programs which are available in File S4.
File S2 Experimental measurements. This folder contains spreadsheets detailing all experimental measurements used in this manuscript. The "rgv2" folder corresponds to CcaSR, "rdv2" to Cph8-OmpR, and "mux" to the dual-system. Within these folders are the spreadsheets, identified using the "Run ID" values indicated in File S1,7-8. Within each of these spreadsheets are a series of tabs corresponding to the experimental measurements (including cytometry histograms, final OD600 values, and incubator temperature time-courses) and the light programs (both in intensity units and in LPA-readable 14-bit grayscale) used to generate them. Each sample can be correlated between tabs using the "Sample ID" column, which provides a unique identifier for each culture sample.

File S3
Reference LED spectra. The raw spectroradiometric reference measurement of each LED is available as a .IRR file (human-readable in a text editor). The post-processed .xls file contains a truncated spectrum and statistics about the LED spectrum. The "reference_leds.xls" file contains the summary statistics for each of the LEDs.
File S4 LED calibration measurements and results. The "process_led_archives.py" Python 2.7 script parses the "led_archives.xlsx" and processes the contents of the "raw" and "sd" folders within each of the LED archive folders according to the LED layouts in the "layouts" folder. The "raw" folder contains the raw calibration measurements, while the "sd" folder contains the SD card files from the LPA used during calibration. The "processed" folder contains spreadsheets with spectra and LED statistics generated for each LED. During processing, an additional file is generated alongside the "processed" folder containing a summary of the LED calibration information for each set of LEDs.
File S5 Reference gene expression signals. The reference programs and perturbation signals used for light program generation are stored in this file. Each reference program consists of a pair of columns identified by "[ID]_x" and "[ID]_y". The first column contains time values (in minutes) while the values of the second column depend upon whether the signal is a reference or a perturbation (as indicated by the ID). Reference signals (ID="ref") describe normalized gene expression levels (i.e. 0 is the minimal output of a system, and 1 is the maximum), while perturbation signals (ID="pert") describe perturbative light programs as a sequence of intensities (in µmol m -2 s -1 ).
File S6 Generated light programs. The light programs produced using the Light Program Generator algorithm are stored in this file. "LED ID" corresponds to the LED which follows the light program, "Times" are the time points (in minutes) corresponding to step-changes in the light intensity, "Intensities" is a sequence of intensities describing the generated light program, "Pre-inten" is the preconditioning intensity used for the program, "Ref ID" is the reference signal used to generate the program, "Model ID" is the identifier of the model used to generate the program, "Compensated perturbation LED ID" is the LED ID for the perturbing LED (if present), and "Compensated perturbation ref ID" is the perturbation signal used by the perturbing LED (if present).
File S7 Plate-by-plate overview of Cph8-OmpR experiments. File contents follow the same description used for File S1. File S8 Plate-by-plate overview of dual-system programming experiments. File contents follow the same description used for File S1.