A photoconversion model for full spectral programming and multiplexing of optogenetic systems

Optogenetic TCS model

We constructed an in vivo optogenetic TCS model comprised of a “sensing model”, which converts light inputs into a ratio of the photoreceptor populations, and an “output model” which converts the photoreceptor populations into a gene expression signal. The sensing model (4) extends the in vitro two‐state photoconversion model to include terms for production of new ground‐state photoreceptors (S_g) at rate k_S and dilution of both 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 and 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 (4) describes a pure time delay (τ) and Hill function mapping between x = 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) (Castillo‐Hair et al, 2016). is the range of possible k_G values, is the minimum value of k_G, n is the Hill coefficient, and K is the 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 GFP production. GFP is degraded in a first‐order process with rate k_dil (4) and has a minimum concentration and concentration range 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⁻² s⁻¹). For such light sources, we can write where is the output spectrum at 1 μmol m⁻² s⁻¹. To quantify the overlap between n_light and σ_i for a given photoreceptor, we introduce as the photoconversion rate per unit light intensity (min⁻¹ [μmol m⁻² s⁻¹]⁻¹). Then, for a given light source, . That is, k₁ and k₂ take on values proportional to light intensity.

Dynamical and spectral characterization of CcaSR

We designed a set of four simple gene expression characterization experiments to train the optogenetic TCS model for CcaSR (Fig 2A–E, Dataset EV2, and Appendix Methods S1 and S2). First, we quantify activation dynamics by preconditioning E. coli expressing CcaSR in the dark, introducing step increases in green light (centroid wavelength λ_c = 526 nm, (Tables EV1, EV2 and EV3, Dataset EV3, and Appendix Method S3) to different intensities, and measuring sfGFP levels over time by flow cytometry (4, Fig 2B, and Appendix Fig S1). Second, we measure deactivation dynamics by preconditioning the cells in different intensities of green light and measuring the response to step decreases to dark (Fig 2C and Appendix Fig S1). Third, we measure the ground‐state spectral response by exposing the bacteria to 23 LEDs with λ_c spanning 369–958 nm with illumination intensities varying over three orders of magnitude (4, Fig 2D, Appendix Fig S2, Tables EV1, EV2 and EV3, and Appendix Method 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 and Appendix Fig S2).

CcaSR model parameterization

Next, we used nonlinear regression to fit the model to the dynamical and spectral characterization data (4, Table EV4, and Dataset EV2). Specifically, we determined and values for each LED, and (LED‐independent) values of the Hill function parameters, k_dil, k_dr, and τ for the system (Fig 2F and G). While simulations using the resulting best‐fit parameters (Fig 2B–E, Dataset EV4, and Table EV4) recapitulate the known properties of the system (Appendix Fig S3), the value of the Hill parameter K is weakly determined. In particular, alterations in K from the best‐fit value can be compensated for by changes in and (Appendix Fig S4). 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 and Appendix Fig S3).

Spectral validation of the CcaSR photoconversion model

Our parameterization experiments yield and values for each calibration LED (Fig 2G). However, to predict the response of an optogenetic tool to a new light source without additional calibration experiments, knowledge of σ_i is required. To estimate σ_i for CcaSR, we developed a procedure to fit a cubic spline to the previously determined and values for each of the 23 LEDs (4, Fig 3A, Appendix Figs S5 and S6, and Datasets EV2 and EV5). 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 σ_i estimate, we measured for a previously untested set of eight color‐filtered white‐light LEDs designed to have complex spectral characteristics (Tables EV1, EV2 and EV3, Dataset EV3, and Appendix Method 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 (λ_c = 526 nm). These predictions are remarkably accurate for LEDs 1–5 (root‐mean‐square errors (RMSEs) from 0.13 to 0.22, 4), which drive sfGFP to high levels, and 7 and 8, which drive low expression (RMSE = 0.17 and 0.14, respectively), but slightly less so for LED 6 (RMSE = 0.28), which drives sfGFP to an intermediate expression level (Fig 3C). These results demonstrate that we can predict the response of CcaSR to a wide range of previously untested light sources using only spectroradiometric measurements of their and not biological calibration experiments.

Dynamic validation of the CcaSR photoconversion model

We previously developed a “biological function generator” method in which we use a predictive model to computationally optimize light input programs to drive tailor‐made gene expression signals such as linear ramps and sine waves (Olson et al, 2014). This 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 and Dataset EV6). 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. We then used the model to computationally design four light time courses each with different LEDs or LED mixtures to program the bacteria to follow this reference signal (4 and Dataset EV6). “UV mono” utilizes a single UV LED (λ_c = 389 nm) (Fig 4A) to demonstrate control of CcaSR with an atypical light source. “Green mono” uses the λ_c = 526 nm LED (Fig 4B) to demonstrate predictive control with a typical light source. “Red perturbation” combines “Green mono” with a strong red (λ_c = 657 nm) sinusoidal signal (Fig 4C and Dataset EV6) designed to demonstrate the perturbative effects of additional sources of light during experiments. Finally, in “Red compensation”, the “Green mono” time course is re‐optimized to compensate for the impact of “Red perturbation” (Fig 4D and 4).

The model predicts the response of CcaSR to all four light signals with high quantitative accuracy (Fig 4 and Dataset EV2). “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 during this experiment (Appendix Figs S7 and S8), and the prediction remains accurate until UV reaches maximum intensity (20 μmol m⁻² s⁻¹). “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 and 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 (Figs EV2–EV4-EV2–EV4, Appendix Figs S9, S10 and S11, Table EV5, and Dataset EV7). 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 in the TCS signaling portion of the output model (4). The model again fits exceptionally well to the experimental data (Fig EV2 and Appendix Figs S9, S10, S11). Unlike CcaSR, which exhibited no detectable dark reversion (Fig 2F), Cph8‐OmpR appears to revert in (Fig EV2F). As before, K is underdetermined (Appendix Fig S4), and we chose the best‐fit value (Table EV5). The Cph8‐OmpR model performs similarly to its CcaSR counterpart in the spectral validation experiments (Fig EV3) and demonstrates greater predictive control in the dynamic validation experiments (Fig EV4).

Click here to expand this figure.

Click here to expand this figure.

Click here to expand this figure.

Development of a CcaSR, Cph8‐OmpR dual‐system model

We engineered a three‐plasmid system (Fig EV1 and Dataset EV1) to express CcaSR and Cph8‐OmpR in the same cell with sfGFP and mCherry outputs, respectively (Fig 5A). 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 sfGFP and mCherry at different combinations of green (λ_c = 526) and red (λ_c = 657) light (Fig 5B, Appendix Fig S12, and Dataset EV8) and refit the Hill function parameters of the TCS signaling portion of the output model (Table EV6). Because the photoconversion parameters are properties of the photoreceptors themselves, we left them unchanged. 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 and Dataset EV8). 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 (4). Due to the minimal response of dual‐system 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 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 EV4), 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 (Datasets EV6 and EV8). “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.

Bacterial strains

All systems utilize the E. coli BW29655 host strain (Zhou et al, 2003). The CcaSR system strain carries the pSR43.6 and pSR58.6 plasmids, which confer spectinomycin and chloramphenicol resistance, respectively (Schmidl et al, 2014). The Cph8‐OmpR system strain carries the pSR33.4 (spectinomycin) and pSR59.4 (ampicillin) plasmids (Schmidl et al, 2014). The dual‐system strain carries pSR58.6, pSR78 (spectinomycin), and pSR83 (ampicillin). Plasmid maps and sequences are available (Fig EV1 and Dataset EV1).

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 (Appendix Method 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₆₀₀ = 0.1–0.2) grown in the same media conditions at −80°C (Appendix Method S2). Cultures were inoculated at low densities (typically OD₆₀₀ = 1 × 10⁻⁵) to ensure that final densities did not reach stationary phase (OD₆₀₀ < 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 LPA instruments (Gerhardt et al, 2016), enabling precise control of two LEDs to define the optical environment of 192 cultures at a time. LPA program files were generated using Iris (Gerhardt et al, 2016) and a custom Python tool (Dataset EV9).

LED measurement

All LEDs were measured and calibrated (Appendix Method S3 and Dataset EV3) 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⁻¹). 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 P_light (λ) (μmol s⁻¹ nm⁻¹) 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 EV1).

Calculation of n_light

Because the LEDs we utilize have fixed spectral characteristics, the spectral flux density (μmol m⁻² s⁻¹ nm⁻¹) incident on the photoreceptors can be parameterized by the LED intensity (μmol m⁻² s⁻¹). 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, n_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⁻¹) divided by the cross‐sectional area of the well. Given a well radius of 7.5 mm, we calculate

LED calibration

Each of the approximately 700 individual LEDs used in the study were measured (Appendix Method S3 and Dataset EV3), enabling compensation for variation in LED and LPA manufacturing (Tables EV1, EV2 and EV3 and Dataset EV9). 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, . 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⁻² s⁻¹ 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⁻² s⁻¹ were re‐measured 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–4,095) 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 (Appendix Method S1) after precisely 8‐h growth by placing the 24‐well 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 (Olson et al, 2014). 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 (Olson et al, 2014). 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 (Castillo‐Hair et al, 2016) via a custom Python script (Dataset EV9). 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 (Castillo‐Hair et al, 2016). 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 previously introduced (Introduction, Results) with best‐fit values summarized in Figs 2F and EV2F. Note that k₁ and k₂ 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 can be expressed as:

where

This ODE can be solved analytically for a step change in light from one environment to another. If the step change occurs at time t = 0, then k₁, k₂, and k_tot are all fixed for t > 0. Given an initial sensor fraction y(0) = y₀, we find.

This solution represents an exponential transition from an initial sensor fraction of y₀ to a final fraction given by with a time constant set by k_tot. As a result, we anticipate that the transition dynamics of y(t) will be slowest under zero illumination when k_tot = k_dil + k_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: k_i = k_{i,source 1} + k_{i,source 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 k_G(t) to the active ratio of light sensors . We model TCS signaling as a pure time delay τ and a sigmoidal Hill function. For CcaSR, the Hill function is activated by increasing sensor ratios, while for Cph8‐OmpR, the inverted TCS signaling activity results in a repressing Hill function. Thus, we write for CcaSR and for Cph8‐OmpR, where the variables and rates have been previously introduced (Introduction, Results) with best‐fit values summarized in Figs 2F and EV2F.

Output gene expression model

We model output gene expression by first‐order production and dilution dynamics:

where the variables and rates have been previously introduced (Introduction, Results) with best‐fit values summarized in Figs 2F and EV2F.

Generation of model simulations

Simulations (Datasets EV2, EV7 and EV8) were produced by numerically integrating the system of ODEs using Python's scipy.integrate.ode method using the “zvode” integrator with a maximum of 3,000 steps (Dataset EV4).

Model parameterization

The CcaSR and Cph8‐OmpR models were parameterized using global fits of the model parameters to the complete training datasets (Figs 2B–E and EV2B–E, and Datasets EV2 and EV7). 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 (Newville et al, 2014). 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 across the complete set of data points .

Estimation of PCSs

Photoconversion cross section estimates were constructed by linearly regressing a cubic spline to the experimentally determined photoconversion rates in order to produce a continuous PCS (Appendix Fig S5 and Dataset EV5). The were produced by minimizing the error between unit experimental photoconversion rates (Figs 2F and EV2F) and spline‐derived predictions . 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 over‐parameterization of the , we used “Leave‐one‐out cross‐validation (LOOCV)” to evaluate the performance of splines with between 5 and 20 knots to determine the ideal number required for each The (Appendix Fig S6). The resulting optimal number of splines was 12 and 8 for CcaS σ_g and σ_a and 12 and 12 for Cph8 σ_g and σ_a. One knot was fixed at 1,050 nm, and the remaining knots were evenly distributed between 350 and 800 nm (Appendix Figs S5 and S6).

Calculation of prediction error (RMSEs)

For model validation, we use a relative error metric ( that reports the root‐mean‐square (RMS) of the log₁₀ error between the predicted and measured responses (Datasets EV2, EV7 and EV8).

Light program generator (LPG) algorithm

The LPG was used as previously described (Olson et al, 2014). The only modification was to use simulations generated by the model described herein rather than the previous model. Compensated light programs were generated by incorporating the presence of the external light signal into the model simulations.

Comparison of output gene expression ranges for single‐ vs. dual‐systems

The CcaSR output range is nearly conserved (60‐fold vs. 56‐fold), while the mCherry response from Cph8‐OmpR is substantially reduced (210‐fold vs. 6.0‐fold). Additionally, the light response is less sensitive than was observed for Cph8‐OmpR individually, as half‐repression requires a 5.2‐fold higher intensity (Appendix Fig S14). We speculate that the reduction in the output range and decrease in sensitivity of Cph8‐OmpR results from a competition between Cph8 and CcaS for limiting PCB, leading to a substantial population of light‐insensitive apo‐Cph8. Notably, the growth rate (Appendix Fig S13) of the dual‐system strain (39.2 min per doubling) is only marginally slower than the single‐system strains (37.4 and 37.9 min for CcaSR and Cph8‐OmpR, respectively).

Detailed descriptions of multiplexed function generation reference signals

In the below descriptions of the multiplexed function generation reference signals (Dataset EV8), the percentages and fractions correspond to a log‐scaled representation of the output range (e.g., if a system has a 16‐fold output range, the 50% level on a log scale would be at the same expression at the 25% level on a linear scale).

1. Dual‐sines. The mCherry reference signal is described by the function while the sfGFP reference signal follows .
2. Sine and steps. The mCherry reference signal is the same as in “Dual‐sines”, while the sfGFP signal is a series of 80‐min holds and 40‐min linear ramps in increasing increments of 20% of the output range.
3. The sfGFP signal is the same as in the “Sine and steps” program, while the mCherry signal is the inverse of the same program.
4. (Time‐shifted waveform). The mCherry signal is a complex function consisting of the following steps:
   1. linear ramp from 0 to 70% over 80 min,
   2. hold at 70% for 40 min,
   3. linear ramp down to 50% over 60 min,
   4. hold at 50% for 40 min,
   5. sinusoidal signal described by the function
   6. hold at 50% for 40 min.

The sfGFP signal is the same program but delayed by 60 min.