# Understanding Cellular Noise with Optical Perturbation and Deep Learning

Chuanbo Liu  State Key Laboratory of Electroanalytical Chemistry, Changchun
Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun, Jilin,
130022, P.R. China.  Yu Fu  State Key Laboratory of Electroanalytical
Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of
Sciences, Changchun, Jilin, 130022, P.R. China.  University of Chinese Academy
of Sciences, Beijing, 100000, P.R. China.  Lu Lin  State Key Laboratory of
Electroanalytical Chemistry, Changchun Institute of Applied Chemistry, Chinese
Academy of Sciences, Changchun, Jilin, 130022, P.R. China.  University of
Chinese Academy of Sciences, Beijing, 100000, P.R. China.  Elliot L. Elson
Department of Biochemistry and Molecular Biophysics, Washington University in
St. Louis, Saint Louis, MO 63130, United States<EMAIL_ADDRESS>Jin Wang
Department of Chemistry and of Physics and Astronomy, State University of New
York, Stony Brook, New York, 11794-3400, United States.
<EMAIL_ADDRESS>
###### Abstract

Noise plays a crucial role in the regulation of cellular and organismal
function and behavior. Exploring noise’s impact is key to understanding
fundamental biological processes, such as gene expression, signal
transduction, and the mechanisms of development and evolution. Currently, a
comprehensive method to quantify dynamical behavior of cellular noise within
these biochemical systems is lacking. In this study, we introduce an
optically-controlled perturbation system utilizing the light-sensitive
Phytochrome B (PhyB) from Arabidopsis thaliana, which enables precise noise
modulation with high spatial-temporal resolution. Our system exhibits
exceptional sensitivity to light, reacting consistently to pulsed light
signals, distinguishing it from other photoreceptor-based promoter systems
that respond to a single light wavelength. To characterize our system, we
developed a stochastic model for phytochromes that accounts for
photoactivation/deactivation, thermal reversion, and the dynamics of the
light-activated gene promoter system. To precisely control our system, we
determined the rate constants for this model using an omniscient deep neural
network that can directly map rate constant combinations to time-dependent
state joint distributions. By adjusting the activation rates through light
intensity and degradation rates via N-terminal mutagenesis, we illustrate that
out optical-controlled perturbation can effectively modulate molecular
expression level as well as noise. Our results highlight the potential of
employing an optically-controlled gene perturbation system as a noise-
controlled stimulus source. This approach, when combined with the analytical
capabilities of a sophisticated deep neural network, enables the accurate
estimation of rate constants from observational data in a broad range of
biochemical reaction networks.

## Introduction

In fully differentiated cells under constant environmental conditions, the
morphology, function, and number of cell-specific organelles, as well as the
capacity for cell-to-cell communication, the potential for interaction with
extracellular molecules, and the function-specific metabolic adaptations,
among other characteristics, appear unchanged from a macroscopic perspective,
suggesting a stable cell state. On the other hand, cell state diversity is
frequently observed through single-cell and single-molecule techniques 1.
Given that the features of a cell are defined by its molecular composition, it
is reasonable to postulate that the number of molecules of proteins, including
enzymes, structural proteins, and signaling molecules, fluctuate around a
constant mean in a stable cell state. These fluctuations are attributed to
both intrinsic and extrinsic sources of noise 2. It is generally assumed that
the magnitude of noise remains consistent over time in a constant environment,
leading to a steady-state distribution. Furthermore, the distinct difference
of living cell to a thermal equilibrium system is the exchange of both energy
and materials with the environment at a constant rate apart from stochastic
fluctuations to maintain the steady-state. From a thermodynamic perspective, a
living cell in this circumstance is characterized as a non-equilibrium steady-
state (NESS) 3.

The ability to reach and maintain NESS hinges critically on the chemical and
thermal equilibrium of the biochemical reactions within the cell, which
defines the phenotype of the cell. Conversely, a cell’s response to external
stimuli is contingent upon the dynamics of signaling components, with encodes
vital biological information for the cell 4, 5. Consequently, the cellular
response to external stimuli is regulated not solely by the steady-state
concentrations of various species but also by the dynamics, which are
characterized by both the instantaneous numbers of molecules and their spatial
distributions. Thus, to thoroughly understand the behavior of biochemical
systems, it is imperative to quantify not just the mean populations of
molecules but, more critically, the temporal variations and spatial
heterogeneity in the distributions of cellular states.

Although significant progress has been made in understanding the noise of
cells in slowly-changing environments 6, 7, 8, 9, the dynamics induced by
rapidly-evolving or spatially-specific perturbations remain poorly understood.
One of the primary challenges is the development of efficient labeling
techniques for molecular species that enable high-throughput tracking and
precise quantification. A further, and perhaps more critical, challenge is the
absence of high-resolution spatio-temporal modulation methods capable of
generating or simulating rapid spatio-temporal stimuli for cellular reaction
networks. The first challenge has been largely addressed through the use of
efficient single-molecule hybridization/tagging techniques or deep sequencing
10, 1, 11. While the second challenge may be surmountable by harnessing
optogenetic technologies. Over recent decades, optogenetics has emerged as a
pivotal technique for exerting precise spatio-temporal control in cells, with
applications in neuroscience, cell physiology, and cardiology 12, 13, 14, 15,
16. Employing a synergistic approach of genetic engineering and bioorthogonal
light, which is independent of the native cellular chemical reactions, both
natural and synthetic photoreceptors can be utilized to modulate the activity
of targeted cells within a highly heterogeneous tissue matrix. The specificity
of expression and subcellular localization within the host organism is
facilitated through the incorporation of regulatory elements such as
promoters, enhancers, and precise targeting sequences into the photoreceptor-
encoding DNA constructs. Consequently, cells that are intrinsically non-
responsive to light can be rendered photosensitive, thereby becoming amenable
to spatially and temporally precise control via targeted illumination.

However, mapping the dynamic properties of biological systems remains
challenging, despite the availability of precise spatial-temporal perturbation
control. From a modeling perspective, constructing a detailed model for gene
expression regulated by light is arduous. The challenge arises because the
rate constants are intertwined with cell-state variables, such as the
populations of RNA polymerase (RNAP) and ribosomes 17, the levels of
transcription factors and micro-RNAs 18, the binding affinities of genes with
RNAP and transcription factors 19, gene copy numbers 20, the phase of the cell
cycle 21, nutrient density 22, and chromosomal conformations 23, among others.
These latent variables, typically unmeasured in biochemical experiments, are
considered intrinsic noise, in contrast to extrinsic noise, which is induced
by the influence of changing environments. The interaction between the control
protocol and cell states, modulated by these latent variables, adds further
complexity to the modeling of the system.

On the other hand, gene expression is an intricate complex chemical reaction
system involving the binding of transcription factors to DNA regulatory
sequences, recruitment of RNA polymerase, synthesis of nascent mRNAs,
modification of these mRNAs, transport of mRNAs from the nucleus to the
cytoplasm, and ultimately, the translation of mRNAs in ribosomes 24. Each of
these processes comprises numerous elementary steps, potentially ranging into
the tens or hundreds 25, 26, 27, 28. Many of these elementary steps remain
poorly understood, either due to the limitations of current technology
rendering them invisible or because they occur at extremely rapid rates 29,
30. Furthermore, acquiring accurate reaction rates for these elementary steps
is particularly challenging. Owing to these complexities, coarse-grained
models have been proposed to simulate the gene expression process. Such models
employ thermodynamic principles, Boolean logic, and differential equations to
deterministically or stochastically represent gene expression, each with
varying degrees of computational efficiency and precision 31, 32. In these
models, gene expression is commonly conceptualized as either a one-state or
two-state chemical reaction system. Transitions between states, as well as
transcription and translation, are often abstracted into effective complex
reactions, each encompassing multiple elementary steps 33. While ordinary
differential equations can provide reasonable estimates of mean population
levels, a stochastic framework offers a more comprehensive approach. This
methodology not only yields insights into the statistical properties of the
system but also generates the population distribution, which is crucial for
understanding the variability inherent in gene expression systems 34, 35.

Furthermore, since stochastic modeling is indispensable, another emerging
challenge is the inferencing of the rate constants, especially for systems
with many species where the state space expands exponentially. This vast state
space precludes the possibility of ergodic integration of the system. To
address this challenge, Monte Carlo methods are often utilized due to their
computational complexity, which is not contingent on the dimensionality of the
state space 36, 37, 38. Despite this, the computational burden remains
substantial. The timescale for an exact one-step simulation of biochemical
reaction networks is typically on the order of milliseconds, whereas
meaningful phenotypic changes occur over hours or days. Consequently,
simulating a single trajectory requires the accumulation of millions of Monte
Carlo steps. In addition, constructing the posterior distribution necessitates
sampling millions to trillions of trajectories, each corresponding to
different sets of rate constants. Even with advanced algorithms and hardware
that facilitate efficient parallel Monte Carlo sampling, inferring rate
constants remains a formidable task, particularly in biological systems with a
large number of species.

In this study, we developed an optically-controlled gene expression system
utilizing the light-sensitive Phytochrome B (PhyB) from Arabidopsis thaliana,
engineered for enhanced spatial-temporal precision to probe the intrinsic
noise within cellular environments. Our system exhibits a continuous response
to pulsed light signals, distinguishing it from previously described blue
light-regulated gene expression systems 39, 40. To broaden the applicability
of our system across various biochemical reaction networks, we constructed a
stochastic model for phytochromes and an integrated model for the entire gene
expression system. These models incorporate the processes of photoactivation,
deactivation, thermal reversion, and the dynamic behavior of the light-
activated gene promoter system. For meticulous control of the system, we
determined the rate constants of the integrated stochastic model using an
omniscient deep neural network. This network is capable of mapping the
combinations of rate constants directly to time-dependent joint state
distributions. By modulating activation rates through alterations in light
intensity and adjusting degradation rates via N-terminal mutagenesis, we
demonstrate that the optically-controlled perturbations can precisely modulate
molecular noise. Our findings highlight the potential of employing an
optically-controlled gene perturbation system as a noise-managed stimulus
source. Furthermore, they suggest the feasibility of inferring rate constants
by monitoring system behavior, facilitated by the predictive capabilities of
the deep neural network.

## Results

### Perturbing gene expression dynamics with light

In order to achieve high spatial-temporal precision in tuning gene expression,
coupled with quantitative control, we have designed an optically-controlled
perturbation system. This system is based on the light-sensitive Phytochrome B
(PhyB) from Arabidopsis thaliana, drawing inspiration from the work of
Shimizu-Sato et al. 41. The full-length PhyB photosensory module (PSM) is
comprised of four domains: the N-terminal extension (NTE), the
Period/ARNT/Single-minded (PAS) domain (nPAS), the cyclic GMP-regulated cyclic
nucleotide phosphodiesterase/Adenylate cyclase/FhlA (GAF) domain that binds
the tetrapyrrole chromophore (P$\Phi$B), and the phytochrome-specific domain
(PHY). At the interface between nPAS and GAF, there exists an intricate
figure-eight knot, referred to as the “light-sensing knot,” which is
intimately associated with the chromophore and plays a critical role in signal
transduction. The PSM’s structure is reinforced by a light-sensing knot loop
(KL) and a hairpin (HP), which interconnect the NTE, GAF, and PHY domains. In
our study, we utilized only the N-terminal fragment (residues 1-621) of PhyB,
encompassing the essential photosensory module (PSM) while lacking the
dimerization/nuclear import module (DNM; also known as the output module in
its prokaryotic counterparts), which facilitates dimerization and nuclear
translocation of the Pfr form 42. In contrast to the full-length PhyB in
planta, the PhyB N-terminal (PhyB(NT)) is monomeric in vitro, devoid of the
DNM 43, 44. For the construction of a light-sensitive, switchable molecular
system, we selected PIF3 as the binding partner for PhyB(NT). PIF3 comprises
three domains: the Active Phytochrome B-binding (APB) motif, the Active
phytochrome A-binding (APA) motif, and the basic helix-loop-helix (bHLH)
transcription factors. The APB motif preferentially interacts with the light-
sensing knot within the N-terminal photosensory module of photoactivated PhyB.
While PIF3 is capable of forming homodimers via its bHLH domain in its full-
length state 45, its interaction with monomeric PhyB(NT) promotes a monomeric
state in PIF3.

For the additional components of the optically-controlled perturbation system,
instead of employing the Gal4 DNA-binding domain (Gal4-DBD) and the Gal4
activation domain (Gal4-AD) fusions for the functional modules, we utilized
the LexA DNA-binding domain (LexADBD) and the VP16 activation domain (VP16AD).
These combinations are commonly employed in yeast two-hybrid assays. We
engineered two fusion proteins, LexADBD-PhyB(NT) and PIF3-VP16AD, to regulate
gene expression via light induction. Both components were integrated into the
yeast genome and expressed under the control of the constitutive promoter
PADH1. After translation, LexADBD-PhyB(NT) is synthesized in the inactive Pr
state and preferentially accumulates in the cytosol in the absence of light.
Consequently, in dark conditions, the photoreceptor accumulates in this
inactive state, ready to be activated upon light exposure. For the reporter
gene construction, LexA binding sites were fused to a segment of the GAL1
promoter encompassing the TATA box. Thus, the PhyB(NT)/PIF3 interaction
modulates expression from a synthetic promoter comprising four LexA binding
sequences, with each sequence capable of binding two LexA molecules. Upon
light activation, the PSM domain of PhyB(NT) undergoes a conformational change
and binds to PIF3-VP16AD. This interaction facilitates the binding of VP16AD
to the TATA box, triggering the expression of a fluorescent reporter for
observation. The degree of photoconversion from Pr to Pfr correlates directly
with the number of red photons administered, which is expected to provide
precise dosage control in a highly predictable and reproducible manner 46.
Light also enables high spatial-temporal precision in control. Such precise
regulation of gene expression is crucial for perturbing cell state dynamics.
To modulate the stability of the fluorescent reporter protein for rapid
degradation and enhanced signal sensitivity, we employed a ubiquitin fusion
technique. This method allows for the alteration of protein degradation rates
by revealing a destabilizing residue following the co-translational removal of
ubiquitin 47. A schematic representation of the optically-controlled gene
expression system is depicted in Fig. 1(A).

The optically-controlled system is expected to respond rapidly in the presence
of red light and to be deactivated by far-red light 48. Owing to the low basal
transcription rate, minimal activation is anticipated in the absence of light
41. We monitored the dynamics of reporter expression using a microfluidic chip
designed for the long-term incubation of thousands of cells (see Methods for a
detailed description of the microfluidic chip fabrication). A representative
experimental result is depicted in Fig. 1 (B). We observed that pulsed light
activation induces sustained protein expression, as illustrated in Fig. 1 (C).
This behavior contrasts with the blue light-regulated gene expression system,
which demonstrated consistent protein expression under continuous light
activation and rapid deactivation in the absence of light 49, 40. Moreover,
this system exhibits significant fluctuations in gene expression compared to
chemically induced gene expression systems. These observations suggest that
protein expression in our system can be finely tuned by modulating the
reaction rates of the light-responsive components and that the underlying
mechanisms of light response in our optical perturbation system are
fundamentally different from those in the blue light-regulated gene expression
system. Prompted by these findings, we proceeded to investigate the
photoconversion and thermal reversion processes of PhyB in the presence of
PIFs.

### Modeling the photoconversion and thermal reversion process of phytochromes

Most land plants, particularly Arabidopsis, possess several phytochromes,
which are prevalent in angiosperms and can be categorized into three primary
types: Phytochrome A (PhyA), PhyB, and PhyC 50. These phytochromes are bilin-
binding photosensory receptors that undergo light-induced conformational
changes. These changes facilitate interactions with signaling partners,
especially transcription factors or the proteins that regulate them, leading
to extensive transcriptional reprogramming. Consequently, phytochromes
orchestrate a variety of physiological processes, including seed germination,
chloroplast development, seedling de-etiolation, shade avoidance,
photoperiodic regulation of flowering time, and ultimately, senescence 51.
Among these types, PhyA and PhyB have similar absorption spectra, yet they
exhibit significant differences in their action spectra for physiological
responses. PhyB is primarily active under red light, while PhyA functions
effectively under far-red light, where PhyB’s activity is suppressed. This
dichotomy has spurred extensive research into their distinct roles.

In the case of PhyB, bilin phytochromobilin (P$\Phi$B) serves as the
chromophore, absorbing energy from red light and transferring it to the
protein’s structure. During the photoconversion process, light-induced bilin
isomerization prompts a conformational rearrangement, converting the dark-
adapted, biologically inactive Pr state into a metastable, far-red light-
absorbing Pfr state, which is biologically active. The partial overlap of
their absorption spectra leads to a dynamic photo-equilibrium between the Pfr
and Pr forms, which varies with the light wavelength. Pfr is produced
exclusively through photoconversion and can revert to the Pr form either
through the absorption of far-red light or via spontaneous relaxation—a
temperature-dependent, non-photochemical process known as thermal reversion.
Thus, the interplay between photoconversion and dark reversion dictates the
steady-state levels of the active Pfr conformation, enabling the dynamic
sensing of light quality and quantity. Consequently, signal transduction is
modulated by a combination of light intensity, wavelength, and temperature, in
accordance with plant requirements 52. Upon exposure to inductive wavelengths,
PhyB adopts a functional dimeric conformation, enabling photosensitive
responses and resulting in three potential phytochrome species: Pr-Pr, Pfr-Pr,
and Pfr-Pfr 53, 44, 42. The intracellular protein dynamics among these PhyB
dimer species dictate their activity potential, with only the nuclear Pfr-Pfr
homodimer initiating PhyB responses 44.

Analysis of in vitro photoconversion induced by red light in Arabidopsis
revealed a monomeric process following single-exponential kinetics 54. In
contrast, the Pfr to Pr photoconversion under far-red light illumination is
more accurately described by a bi-exponential kinetic fit. Additionally, at
room temperature (22∘C), total Pfr reverted entirely to the inactive form (Pr)
following bi-exponential decay kinetics in darkness 43, 54. The bi-exponential
kinetics observed in both photoconversion and thermal reversion may result
from sequential reactions, the presence of multiple conformers, or a
combination of these scenarios. Previous studies have identified intermediates
of phytochromes on the millisecond timescale 55, 56, yet no intermediates have
been detected on the timescale of hours during the thermal reversion of
phytochromes. Consequently, we propose that the bi-exponential decay arises
from two distinct forms of Pfr within the molecular population 57, 58, 59, 60,
43, 52, 42. On the other hand, thermal reversion is also modulated by protein
interactions. The Pfr-Pfr homodimers exhibit greater stability than the Pfr-Pr
heterodimers; the additional subunit interactions in Pfr-Pfr homodimers confer
a stabilization of approximately 4 kcal/mol 61. Photoconversion of only one
subunit in the phytochrome dimer from Pr to the metastable Pfr conformation
could lead to a reduced stability of the Pfr-Pr heterodimer, thereby favoring
a light-independent back reaction to the Pr-Pr homodimer. This phenomenon
results in a more pronounced decline of PhyB responses at longer wavelengths.
The rate of thermal reversion for the Pfr-Pr heterodimer is approximately
100-fold faster than that of the active Pfr-Pfr homodimer, enabling PhyB to
perceive temperature fluctuations both diurnally and nocturnally 44. Studies
of Arabidopsis PhyB expressed in yeast have revealed that the dark reversion
rates of Pfr-Pfr homodimers are an order of magnitude lower than those of Pfr-
Pr heterodimers 61. The interaction of PhyB with ARR4 (Arabidopsis Response
Regulator 4) diminishes the dark reversion of the Pfr conformer, thus
modulating light sensitivity in Arabidopsis 62. Furthermore, phytochromes
demonstrate prolonged persistence of the Pfr state when associated with
nuclear bodies (NBs), attributable to the suppression of thermal reversion in
darkness 44. The binding of PhyB(NT) to PIF can also inhibit the thermal
reversion process 43, 52. Collectively, these findings suggest that the
photoconversion rate of phytochrome from Pfr to Pr is contingent upon the
wavelength and fluence rate of light, while the thermal reversion rate is
influenced by temperature and the presence of binding partners.

To consolidate the findings of prior studies, we introduce a streamlined model
to encapsulate the photoconversion and thermal reversion mechanisms of
phytochromes, as depicted in Fig. 2 (A). This model is akin to the dual-
pathway framework suggested by 63. It posits the existence of two photo-
activated conformers, Pfr1 and Pfr2, inferred from the bi-exponential kinetics
characteristic of the thermal reversion process 54. The model allows for the
PIF binding partner to interact with both the Pfr and Pr states, with the
latter also capable of basal affinity binding to Pa 64. Photoconversions are
postulated to occur in the dimeric state, with rate constants mirroring those
of the monomeric state. To determine these rate constants, we devised a score-
based bootstrap inference algorithm that leverages in vitro equilibrium data
from photoconversion and thermal reversion assays, as illustrated in Fig. 2
(B, C). The score is defined by the equation

$\mathcal{S}=\sum_{t}\int\left\|[\boldsymbol{x}(t)]-\boldsymbol{y}(t)\right\|_{2}d[\boldsymbol{x}(t)],$
(1)

representing the mean L2-norm discrepancy between the model’s theoretical
predictions and the empirical data across all observable time points. For each
set of experimental conditions, we compute the score using 10 simulated
trajectories, each generated with the Gillespie exact kinetic sampling
algorithm 37. The inference algorithm iteratively adjusts the rate constants
to minimize the score via the steepest descent approach. The photoconversion
rates are influenced by the light’s wavelength and fluence rate; red light
catalyzes both the Pr $\to$ Pfr and Pfr $\to$ Pr transitions in phytochromes.
For the sake of simplicity, we characterize photoconversion using effective
rates corresponding to specific light wavelengths. That is, when modeling the
Pr $\to$ Pfr (Pfr $\to$ Pr) photoconversion, the rate for the Pfr $\to$ Pr (Pr
$\to$ Pfr) transition is set to zero. Additionally, while Arabidopsis PhyB(NT)
exhibits single-exponential kinetics under red light, we assign different
photoconversion rates to Pr $\to$ Pfr1 and Pr $\to$ Pfr2 to achieve a more
accurate model, as evidenced by a reduced score.

The chemical reactions and rate constants inferred are detailed in the SI
Appendix, Table S2. The photoconversion of Pfr1 $\to$ Pr occurs approximately
3-fold faster than that of Pfr2 $\to$ Pr, while the thermal reversion of Pfr2
$\to$ Pr is roughly two orders of magnitude slower than that of Pfr1 $\to$ Pr.
Consequently, the Pfr2 state is a photo-activated state that forms slowly but
exhibits greater stability compared to Pfr1. Upon binding to PIFs, the thermal
reversion rate is reduced by an order of magnitude for both Pfr1 $\to$ Pr and
Pfr2 $\to$ Pr, with the reversion of Pfr2 $\to$ Pr remaining two orders of
magnitude slower than that of Pfr1 $\to$ Pr. The dissociation constants
$K_{D}$ for Pfr1Pa and Pfr2Pa are 0.56 and 23.81, respectively, signifying
that Pfr1 has a stronger affinity for PIFs than Pfr2. For PrPa,
$K_{d}=6.23\times 10^{4}$, which is approximately four orders of magnitude
greater than that for Pfr1 and two orders of magnitude greater than that for
Pfr2, aligning with the findings of previous studies 64. Our analysis suggests
a model with two photo-activated conformers that can account for the observed
photoconversion and thermal reversion behaviors of phytochromes both in
isolation and when associated with PIFs. The photoconversion rates of Pr $\to$
Pfr1 and Pr $\to$ Pfr2 are comparable. However, Pfr2 demonstrates enhanced
stability in terms of both photoconversion and thermal reversion. The binding
of PIFs further stabilizes the photo-activated state, leading to a reduced
rate of deactivation that mirrors the behavior observed in planta 61, 62, 44,
65, 66, 42. As depicted in Fig. 2 (C), our findings surpass those of the prior
study by Smith et al., which did not accurately predict the thermal reversion
process at elevated PIF concentrations 43. Moreover, our model presupposes the
existence of only two activated conformers, obviating the need to postulate a
secondary activated state for phytochromes. Collectively, these results affirm
that our model effectively captures the comprehensive behavior of
photoconversion and thermal reversion processes in phytochromes, both
independently and in the presence of PIFs. We have therefore employed the
framework of this streamlined model to analyze the optically controlled gene
expression system.

### Stochastic modeling and parameter screening with omniscient neural network

We monitor the real-time kinetic signal to conclude that the system has a
distinct dynamical behavior compared to those systems activated by blue light
49, 40. In the blue light-activated systems, the light-switchable promoter
gradually ceases to produce protein upon light removal. However, our system
behaves differently: it is activated and remains active until exposed to
deactivating light. The light filtering pattern in the blue light-activated
systems can originate from its light-out dissociation (LOD) behavior. In
contrast to these systems, where the memory of the expression state is
gradually lost in darkness, the dimeric dynamics of PhyB lock the system in
the activated state, leading to a robust memory effect. Considering both
strong and weak memory, and the confidence in these two mechanisms, the
biophysical properties of our system endow it with unique signaling
capabilities: a strong memory signaling with an enhanced transfer function.

Given the highly stochastic fluorescent signals from the gene expression
trajectories controlled by PhyB-based optical regulators (Fig. 1 (C)), a
deterministic description of the system is precluded 34, 35. Therefore, from
the perspective of probabilistic modeling, we represent the optically-
controlled gene expression regulated by phytochromes as a non-Markovian
stochastic dynamical system that is responsive to external controls. A
stochastic reaction network has been constructed for this system 67. In light
of the evidence presented in the previous section, which reveals two activated
conformers through fitting to photoconversion and thermal reversion processes,
we have modeled the system comprising 23 species and 74 reactions, as depicted
in Figure 3(A) and listed in Table S3.

The control of gene expression by light, as described in the previous section,
involves altering light exposure to red or far-red light, which modulates the
rate constants of photoconversions and can initiate or terminate mRNA
transcription. The activation function of the promoter is modeled as a step
function, aligning with the timescale separation where PhyB-PIF3 binding
occurs within milliseconds, and gene transcription takes place over seconds.
To address the fluorescent background issue associated with the basal
expression of the fluorescent reporter, we initiate the modeling 12 hours
prior to the experimental start time. The complete control protocol is as
follows: 12+2 hours of dark incubation, a 2-minute red-light pulse for
activation, a subsequent 5-hour dark incubation, and finally an 8-hour
deactivation period under far-red light. The system then reaches a NESS prior
to the experiment, reflecting the actual system setup. The background cellular
fluorescence signal intensity is then used to calibrate the steady-state
distribution of the fluorescent reporter generated by basal expression.
Consequently, the distribution of the background signal also provides insights
into the rates of basal transcription and degradation. Furthermore,
calibrating the background signal determines the ratio factor $R$ for
converting fluorescent intensities to the number of molecules, defined by the
equation:

$R=\frac{\text{Fluorescent intensity (a.u.)}}{\text{Number of reporter
molecules}}$

where $R$ is expressed in arbitrary units (AU). During the experimental time
frame, the ratio factor is assumed to remain constant, independent of the
number of reporter molecules present. Hence, for each experiment, a consistent
ratio factor is applied to accurately fit the stochastic model to the
experimental data.

For extended two-state stochastic gene expression systems, distributions can
be analytically obtained 33. However, solving the PhyB-based optically-
controlled gene expression system analytically is infeasible due to its large
state space. Numerical methods, such as kinetic Monte Carlo methods 37, can be
employed to approximate the joint state probability distribution. Yet,
sampling from the posterior distribution given a general prior and
experimental data remains impractical with these methods because of the
prohibitive computational time. For example, estimating the likelihood
requires simulating approximately $10^{4}$ trajectories for each set of rate
constants. Typically, for a single set of rate constants, each estimation
takes about 10 minutes using a parallel algorithm. To accurately estimate the
posterior distribution, it is necessary to sample $10^{4}$ combinations of
rate constants, not including the potentially more numerous mixing steps.
Consequently, even with a reasonably informative prior, one inference can take
an estimated 60 days to complete, rendering the process highly inefficient.
This limitation is analogous to other numerical methods, such as finite state
projection approaches, which track the joint state probability distribution
within a truncated state space 68.

To surmount this obstacle, we employ a neural network strategy to estimate the
time-dependent state distribution using variational approximations. The neural
network’s architecture is depicted in the SI Appendix, Fig. S2. The inputs to
this network consist of prompts, which include the rate constants
$\boldsymbol{\sigma}$, the initial condition $\boldsymbol{x}_{0}$, and the
time value $t-t_{0}$, as well as the state vector at the current time
$\boldsymbol{x}_{t}$. The network’s outputs are the conditional marginal
probabilities of the state coordinates given the prompts’ information, denoted
as
$\hat{p}^{\boldsymbol{\theta}}(x_{t}^{(i)}|x_{t}^{(1)},\ldots,x_{t}^{(i-1)},\boldsymbol{\sigma},\boldsymbol{x}_{0},t-t_{0})$.
The conditional joint probability of the state is thus expressed as

$\hat{p}^{\boldsymbol{\theta}}(\boldsymbol{x}_{t}|\boldsymbol{\sigma},\boldsymbol{x}_{0},t-t_{0})=\prod_{i=1}^{N}\hat{p}^{\boldsymbol{\theta}}(x_{t}^{(i)}|x_{t}^{(1)},\ldots,x_{t}^{(i-1)},\boldsymbol{\sigma},\boldsymbol{x}_{0},t-t_{0})$

To delineate the mapping from the amalgamation of prompts and states to the
conditional marginal probabilities, we implement a masked multi-head self-
attention block that computes the scaled dot-product attentions. The outputs
from the multi-layer decoder, composed of masked multi-head attention blocks,
are subsequently fed into a language model head. This head is designed to
project the decoder outputs onto probability weights, which are then
normalized using the softmax function (see Methods for more details).

To demonstrate the accuracy and efficiency of the neural network approach, we
trained a neural network on a dataset comprising 256 trajectory ensembles,
each consisting of 1,000 trajectories corresponding to a unique combination of
rate constants. The trained neural network subsequently provided the
conditional likelihood
$\hat{p}^{\boldsymbol{\theta}}(\boldsymbol{x}_{t}|\boldsymbol{\sigma},\boldsymbol{x}_{0},t-t_{0})$,
facilitating parameter screening. As illustrated in Fig. 3 (B), despite the
neural network being trained on only 256 trajectory ensembles, the means and
standard deviations of states generated from the neural network’s predictions
for 2,400 prompts closely mirror those obtained from stochastic simulations.
Owing to the highly parallel computation capabilities of the neural network on
General-Purpose Graphics Processing Units (GPGPUs), the computation of the
joint state probabilities for 2,400 prompts can be completed within one hour.
This represents a significant improvement in efficiency compared to
traditional methods, which typically require approximately 40 hours for
similar computations.

### Inferencing the rate constants with neural network assisted Monte Carlo
sampling

Estimating the likelihood using a trained neural network provides a unified
framework for inferring rate constants from a dataset that includes both time-
slice and kinetic data (see Methods for more details). This approach allows
for the integration of information derived from multiple measurement
techniques, encompassing both time-slice and real-time kinetics of any
species, to enhance the accuracy of inference. In contrast, the previous non-
parametric Bayesian method utilized only time-slice data for a limited number
of species 38. For systems with many species, solving the chemical master
equation is computationally formidable, and representing the state
distribution is also memory prohibitive. Nevertheless, our neural network-
based method is not constrained by these limitations (unpublished data).

Given the constraints of computational resources, we have approximated certain
rate constants using values sourced from prior researches. The basal
transcription rate of the reporter gene, absence from any regulatory binding,
was established at 0.001 min-1. This rate negligibly impacts the basal
reporter distribution. Transcription rates for PhyB and PIF were inferred from
the nascent transcription rates of yeast genes, as reported in 69, 70, which
vary between 2 and 30 hour-1. Adopting the upper limit of this range, 0.5
min-1, we approximated the transcription rates of PhyB and PIF to be 0.31
min-1 and 0.33 min-1, respectively. We postulate identical activated
transcription rates for the reporter gene when bound to Pfr1Pa, Pfr2Pa, and
PrPa. Translation rates for PhyB, PIF, and the reporter were estimated from
translation initiation and elongation rates derived from ribosome profiling
data 71, with the translational rate constant ranging between 5 and 8 residues
sec-1. By applying the highest observed rate, we estimated the translation
rates for PhyB(NT)-lexDBD, PIF3-VP16, and the reporter inclusive of a
ubiquitin fragment to be 9.573 min-1, 0.775 min-1, and 1.89 min-1,
respectively. Post-translation, proteins must properly fold to become
functional. The median folding time, $t_{50}=0.005$ sec, as sourced from a
protein database and adjusted for temperature, equates to a rate constant of
8317.8 min-1 is used for both PhyB(NT)-lexDBD and PIF3-VP16 72. The maturation
time for the reporter (sfGFP), $t_{50}=6.9$ min, was derived from in vivo
studies of fluorescent protein maturation in yeast cells, corresponding to a
rate constant of 0.1 min-1 73. We have assumed uniform mRNA degradation rates
for PhyB, PIF, and the reporter. Additionally, we posit that the degradation
rates for both immature and mature proteins are equivalent.

Four rate constants were identified as having the most significant effects on
the behavior of the system. The first is the activated transcription rate of
the reporter ($\sigma_{PR}$), which largely determines the steady-state
distribution of the number of molecules. Another is the association rate of Pr
with Pa ($\alpha_{PrPa}$), which greatly influences the basal expression
level. Finally, the degradation rates of mRNAs ($\beta_{mRNA}$) and the
reporter ($\beta_{PR}$) define the steady-state distribution after
deactivation. Using the trained neural network model, we drew samples from the
posterior distribution of the rate constants using a Markov Chain Monte Carlo
method (see Methods for details). As shown in Fig. 3 (C), the four rate
constants, $\sigma_{PR}$, $\alpha_{PrPa}$, $\beta_{mRNA}$, and $\beta_{PR}$,
are tightly constrained by the experimental data. Utilizing the most probable
values of these rate constants for inference, we performed a stochastic
simulation of the system and compared the predicted state distribution with
the experimental histogram, as depicted in Fig. 3 (D). The complete plots
across 48 time points are presented in the SI Appendix, Fig. S5. It is evident
that the inferred rate constants accurately capture the state distribution of
the fluorescent reporter. The complete list of rate constants is provided in
the SI Appendix, Table S3.

### Modulation of gene expression noise by light

For the purpose of applying our optical perturbation system to the study of
cellular noise in biological systems, it is essential to be predictive for the
output noise of our system. Thus, we investigated the possibility of precise
noise control by employing variable light intensities and modulation of the
degradation rate of the target molecule. To tune the exposure light
intensities, we calibrated the light source of the microscope and applied a
linear increase in light intensities to activate the system. Additionally, we
modulated the degradation rate of the reporter protein by expressing a
chimeric ubiquitin-X-reporter construct. Upon expression, ubiquitin is co-
translationally removed from the protein, leaving behind the X-reporter, where
X represents the chosen amino acid residue. The N-terminal residue of a
protein significantly influences its half-life, which can vary from several
hours (as in the case of methionine or alanine) to a few minutes (for residues
like leucine or arginine) 74, 75. Subsequently, we evaluated the behavior of
our system under conditions of linearly increasing light intensities and
following N-terminal mutagenesis.

As shown in Fig. 4 (A, B), the steady-state distributions of the reporter
indicate changes in the rate constants due to varying light intensities and
N-terminal mutagenesis. We further sought the optimal rate constants for
photoconversions as well as degradation using a score-based bootstrap
inference algorithm (Eq. 1). The comparisons between the molecular number
distributions predicted by the inferred rate constants and the experimental
histograms is presented in the SI Appendix, Figs. S3–S7. The Hellinger
distances between the predicted distributions and the experimental results are
around 0.1, which attests to the efficacy of the inference algorithm. From
Fig. 4 (C), it is observed that within the range of 20%–60% light intensities,
there is an almost linear increase in the rate constant for the reaction Pr
$\to$ Pfr2. In contrast, the rate constant for the reaction Pr $\to$ Pfr1
increases exponentially. Consequently, the extent of photoconversion from Pr
state to Pfr states is directly correlated to the quantity of red photons
provided, exhibiting high predictability and reproducibility. However, both
rate constants diminish to minimal values when the light intensity reaches
80%, likely due to phototoxicity at high light fluences 41.

To quantify the change in cellular noise during the modulation of rate
constants, we calculated the logarithm of the coefficient of variation (CV),
defined as $\log_{10}(\sigma/\mu)$, where $\sigma$ represents the standard
deviation and $\mu$ denotes the mean number of molecules. As depicted in Fig.
4 (D), the logarithm of the CV decreases with an increase in activation rates
and increases with the degradation rate. The dependence of the logarithm of
the CV on the rate constant of Pr$\to$Pfr2 is linear, whereas its dependence
on the input light intensities is predominantly exponential. The slight
deviations from linear and exponential dependencies suggest that the Pr $\to$
Pfr2 transition contributes significantly to the noise behavior of the system.
Furthermore, an exponential dependence on light intensities implies that the
noise intensity has an intrinsic lower bound in the PhyB-based optical
perturbation system. Conversely, the noise’s dependence on the degradation
rate constant is exponential, indicating that the degradation rate constant’s
influence on noise is pronounced when its value is small.

The complex behavior of the system can be comprehensively understood through
the lens of the so-called inverse correlation between noise and mean
expression level, as illustrated in the SI Appendix Fig. S8. Both activation
and deactivation in the optically-controlled perturbation system modulate the
expression level and the expression noise simultaneously, with the noise being
linearly dependent on the mean population of molecules. This relationship
between noise and mean population level has also been observed for both
protein and mRNA in numerous other systems 76, 77, 78, 79. However, the PhyB-
PIF3 system exhibits a high level of noise relative to its mean expression
level compared to other light-modulation systems 76. Considering the dual-
pathway mechanism of the PhyB system, it is conceivable that an increase in
transition pathways maintains the linear dependency relationship while
elevating the ratio of this linear dependency between noise and the mean
population level of molecules.

## Discussion and conclusion

Understanding the behavior of noise in biological systems is crucial for
comprehensively characterizing the phenotypic changes in many key processes,
including gene expression, signal transduction, development, and evolution.
While numerous single-molecule labeling methods are available and high
spatial-temporal resolution can be achieved through optogenetic techniques,
quantitatively describing noise remains a significant challenge. Two major
obstacles are the modeling of complex biochemical reactions and the inference
of rate constants for these models, both of which have garnered considerable
attention.

In this work, we present a systematic approach to elucidate the dynamic
behavior of cellular noise using our optical perturbation system combined with
deep learning. We demonstrate that the optically-controllable system we have
developed constitutes a highly sensitive light modulation system. It exhibits
distinct behaviors compared to the blue light-regulated gene expression system
49, 40. The dynamic mechanism of the blue light-regulated system can be
characterized by Hill functions, whereas our PhyB-based system is more
accurately described by a multiple-path reaction system involving two
activated conformers. Moreover, the blue light-regulated system has a faster
thermal reversion rate compared to the PhyB-based system. Activation by blue
light of these systems leads to a gradual deactivation in these systems,
resulting in a decreased population of the activated state. However, the
stability of the activated Pfr states of phytochromes is significantly
enhanced in the presence of PIFs. As the result, the PhyB-based system not
only retains memory for long time duration but also allows for the
accumulation of light fluence over prolonged periods. This characteristic
renders it a highly stable two-state light-regulated gene regulator when
compared with the blue light system. Both pulsatile and continuous activation
modes can be achieved through the interplay of red and far-red light, without
the need for continuous light illumination. This endows the PhyB-based system
with advantages in in vivo applications, such as the semi-automatic control of
glucose homeostasis or immunotherapy by epigenetically engineered
photosensitive cells 80, 14.

To further quantitatively characterize the dynamical properties of our system,
we developed a stochastic model comprising 74 reactions. To infer the rate
constants of this model, we introduced an omniscient neural network capable of
mapping the rate constants directly onto joint state distributions. The
trained network facilitates efficient computation of likelihoods for various
rate constant configurations. We demonstrate that rate constants can be
determined with high efficiency and precision using the omniscient neural
network, surpassing conventional methodologies. Comparing the stochastic
model’s predictions with experimental data, we show that our system can not
only produce a controllable level of target molecules but also modulate
expression noise through adjustments in light intensity and the degradation
rate of the target protein.

In conclusion, we showcase a novel integration of optical perturbation with a
neural network-based likelihood estimator to dissect noise in biological
systems. Our methodology is readily extendable to diverse optogenetic systems
and organisms, enhancing the elucidation of complex dynamical behaviors in
state distributions for expansive systems.

## Material and methods

### Strains and plasmids construction

All experiments were carried out in W303a (MATa leu2-3,112 trp1-1 can1-100
ura3-1 ade2-1 his3-11,15), a haploid derivative of W303 obtained from ATCC
(#208352). Plasmids were constructed and propagated in E. coli DH5$\alpha$ and
verified via Sanger sequencing prior to strain construction. Standard
protocols for the growth, maintenance, and transformation of yeast and
bacteria, as well as DNA manipulation, were employed. Initial plasmids were
assembled as described in a previous study 41, with functional fragments
cloned from cDNA by conventional PCR. DNA digestion and ligation utilized
restriction endonucleases and T4 DNA ligase from New England Biolabs (USA).
Light-responsive strains were generated by integrating pAM206 into the trp1
locus of the W303a yeast genome. Given that pAM206 harbors two XbaI sites,
partial digestions were executed using the XbaI endonuclease to optimize the
yield of linearized plasmids. These linearized plasmids were subsequently
transformed into yeast cells following the standard lithium acetate-based
transformation protocol. Transformed cells were cultured on selective
auxotrophic minimal media using the appropriate dropout plates (Takara).
Integration of plasmids was confirmed by genomic PCR at both the N- and
C-terminal sides of the recombination loci. In light of the fact that multiple
copies can lead to high expression levels and considerable expression noise,
potentially complicating expression kinetics 76, single-copy integration was
verified using real-time quantitative PCR. The resultant strains expressed
PhyB(NT)-lexADBD. Subsequently, pAM205 was linearized with ClaI and integrated
into the leu2 locus, resulting in strains that additionally expressed
PIF3-VP16AD. The strain produced through sequential integration served as the
progenitor for all subsequent experiments. To modify the stability of the
reporter, pAM236M was engineered by appending an ubiquitin sequence upstream
of the fluorescent protein in pAM236. The N-terminal methionine (M) of pAM236M
was altered to arginine (R) through site-directed mutagenesis, thereby
modulating the stability of the reporter via the N-end rule pathway 74. To
replace the fluorescent reporter Venus with sfGFP, which offers enhanced
photostability and accelerated maturation 81, the sfGFP sequence was initially
PCR-amplified from the GTTg plasmid (Addgene #81109) and subsequently fused
with the backbone sequences of pAM234 and pAM236 via inverse PCR.
Comprehensive plasmid information is available in Table S1. The functionality
of the final constructs was confirmed by fluorescence microscopy following co-
incubation with PCB. Control assays were conducted using PhyB-only and
PIF3-only cell strains under identical conditions for validation.

### Fabrication of microfluidic chip for real-time kinetic measurements

The design and fabrication of the microfluidic device for live-cell imaging
adhered to protocols from previously published work 82. Briefly, the process
utilized SU-8 2000 series photoresists (MicroChem), chrome glass masks (HTA
Photomask), and an EVG620 contact mask aligner (EV Group) to construct and
pattern the desired features on silicon wafers supplied by University Wafer
Inc. The feature heights were confirmed with a Dektak 150 surface profiler
(Veeco). Each chip comprised four identical microfluidic devices, each
inspected meticulously to ensure the absence of dust or defects. Prior to cell
introduction, the microfluidic chip was cleansed with acetone and isopropanol,
followed by at least two rinses with sterilized deionized water to remove any
residual solvents. The chip was then positioned on a custom gas chamber
situated on the stage of an inverted microscope. This setup was enclosed
within a thermally stable microscope cage incubator (Okolab), maintaining a
constant temperature of 30 ∘C. The media ports were connected to plastic
tubing, which in turn was linked to a high-precision pump supplying fresh
media. The waste port of the chip was attached to plastic tubing leading to a
50 mL collection tube.

### Media and growth conditions

Cells were initially cultured in synthetic complete (SC) medium (pH 6.0)
overnight until the optical density at 600 nm (OD600) reached approximately
1.0, prior to initiating experiments. For cell loading, the yeast culture was
transferred using a syringe, driven by a pump at a flow rate of 0.1 mL/min.
The majority of the traps were populated with cells within 1-2 minutes.
Subsequently, the medium in the syringe was replaced with SC medium
supplemented with 25 $\mu$M of PCB, and the cells were incubated for a minimum
of 2 hours in the absence of light. For the light-activated kinetic assays,
the flow rate of the medium was set to 0.01 mL/min to minimize cell movement
and maintain a normal growth rate.

### Light activation procedure and measurements

The structural analogs of phytochromobilin (P$\Phi$B), phycocyanobilin (PCB)
(CAS No. : 20298-86-6), are utilized to induce the photoconversion of PhyB. We
employed a Semrock 641/75 nm bandpass filter and a Semrock 715 nm long-pass
filter for light activation and repression, respectively. For light
activation, red light at 20%, 40%, 60%, and 80% constant intensities was
applied for 2-minute intervals to initiate gene expression. Conversely, for
light repression, far-red light at a constant 50% intensity was used to
effectively suppress transcription in the light-regulated system. Fluorescent
reporter activity was monitored using 488 nm laser illumination with an
exposure time of 100 ms.

### Image analysis with point-prompt cell segmentation

To analyze the microscopic images, we have developed a Python-based software
that facilitates efficient cell segmentation and cell tracking. The software
employs the Segment Anything Model (SAM) as its core model, which is capable
of segmenting yeast cells using point-prompt information 83. Following the
SAM-based segmentation, the masks of the segmented cells are manually refined
through interactive labeling. Subsequent to segmentation, cell labels are
reassigned using automatic tracking via the Hungarian algorithm. The
automatically segmented and tracked cells are subsequently manually reviewed
to correct any mislabeling and to exclude dead cells. Finally, each cell’s
fluorescence intensity, as well as morphological information, are collated and
exported to a text file for additional analysis.

### Simulation of stochastic reaction network

To account for the basal expression of reporter proteins, we simulated the
entire process from -12 hours until the initiation of the experiment. Red
light activation was applied at 840 minutes for 2 minutes, and the system was
deactivated using far-red light at 1142 minutes until the conclusion of the
experiment. The total duration of the simulation was 1624 minutes. We utilized
the biocircuits package for parallel simulations employing the Gillespie
algorithm. For each set of rate constants, $10^{3}$ – $10^{4}$ trajectories
were generated, and system states were recorded at time points corresponding
to the experimental sampling times. To initially train the autoregressive
model in a supervised manner, we simulated trajectory ensembles for 256 sets
of rate constants, where the rate constants were sampled using the quasi-Monte
Carlo method with the Sobol sequence.

### Architecture of the omniscient neural network model

For modeling the light-regulated gene expression system, we have used the
variational approach to approximate the joint distribution of states given
prompts that encode information of rate constants and time. We designed the
omniscient neural network model with masked multi-head self-attention blocks
as its core functional module which describes the connection between prompts
and states. For calculating the self-attention score, we use learnable
embeddings to convert the state vector $\boldsymbol{x}\in\mathbb{R}^{N}$ to a
tensor in the embedding space $X\in\mathbb{R}^{N\times d_{emb}}$. On the other
hand, we collect the logarithm of rate parameters and time value in sequence
to construct the prompt. To handle the continuous variables of rate parameters
and times in the prompt, we first transform rate parameters into logarithmic
values to reduce value divergences. Then, we employ a one-layer perceptron to
project the prompt to a vector $\boldsymbol{p}$ with fixed length $d_{p}$. We
then map $\boldsymbol{p}$ to the same embedding space as the discrete-valued
states with linear projector, to obtain the tensor
$Y\in\mathbb{R}^{d_{p}\times d_{emb}}$. The concatenation of $X$ and $Y$
serves as the input tensor $Z=\text{Concat}(X,Y)\in\mathbb{R}^{d_{s}\times
d_{emb}}$ to the network, where $d_{s}=N+d_{p}$. The masked multi-head self-
attention block calculates the scaled dot-product attention of the input
tensor $Z$ with $h$ attention heads

$\displaystyle\text{Attention}(Z)$
$\displaystyle=\text{Concat}(\text{head}_{1},\ldots,\text{head}_{h})W_{O}$
$\displaystyle\text{where}\quad\text{head}_{i}$
$\displaystyle=\text{Softmax}\left(\frac{\text{Mask}(Q_{i}K_{i}^{T})}{\sqrt{d_{k}}}\right)V_{i}$

where the query $Q_{i}=ZW_{i,Q}\in\mathbb{R}^{d_{s}\times d_{k}}$, the key
$K_{i}=ZW_{i,K}\in\mathbb{R}^{d_{s}\times d_{k}}$, and the value
$V=ZW_{i,V}\in\mathbb{R}^{d_{s}\times d_{v}}$ are all linear projections of
the input tensor $Z$. These projections are obtained using different learnable
parameter matrices: $W_{i,Q}\in\mathbb{R}^{d_{emb}\times d_{k}}$,
$W_{i,K}\in\mathbb{R}^{d_{emb}\times d_{k}}$,
$W_{i,V}\in\mathbb{R}^{d_{emb}\times d_{v}}$, and
$W_{O}\in\mathbb{R}^{hd_{v}\times d_{emb}}$. Here, $d_{k}$ and $d_{v}$
represent the dimensions of the number of query (key) elements and value
elements for each state value, respectively. When $h$ is specified, $d_{k}$
and $d_{v}$ are defined according to the relation $d_{k}=d_{v}=d_{emb}/h$. The
attention mask $\text{Mask}(\cdot)$ is used to conceal the dependency
relations between the previous position and subsequent positions of both the
prompt and the generated sample states. This is achieved by assigning
$(Q_{i}K_{i}^{T})_{\xi<\mu}=-\infty$, where $\xi$ and $\mu$ are the row and
column indexes of the $Q_{i}K_{i}^{T}$ matrix, respectively. We stack multiple
attention blocks to create a multi-layer decoder, which enhances network’s
representational capability. The output of the masked multi-head attention
block $A\in\mathbb{R}^{d_{s}\times d_{emb}}$ is used as the input for another
block. Following the multi-layer decoder, a language model head is utilized to
map the decoder outputs to probability weights, which are then normalized
using the softmax function. The predictions of the MET provide joint
distributions of states given the input prompts.

### Coordinate embedding and training of the omniscient neural network model

The synthesis of mRNA and the expression of proteins are significantly
influenced by the biochemical rates associated with each reaction. The number
of protein molecules ranges from hundreds to millions, necessitating a large
state space for neural network representation. To reduce the number of states,
we employed coordinate embedding to transform state coordinate values into
vectorized representations. This embedding process involves base conversion of
the coordinate value into a binary or decimal number, followed by the creation
of a vector with each element corresponding to a digit in the number. In the
case of a light-controlled gene expression system, to achieve efficient
computations, we constrained the vector sizes for certain coordinates to
automatically account for state constraints, such as the number of genes. This
reduction in vector sizes also decreased the number of species. Training was
conducted using a trajectory ensemble dataset comprising 0.25 million
trajectories. To generate this ensemble, 256 different parameter sets were
initially sampled using a quasi-Monte Carlo method with the Sobol sequence.
Subsequently, 1,000 trajectories were sampled for each parameter set using
Gillespie’s algorithm. The model underwent training for 40 epochs, totaling 80
hours on a cloud server equipped with 4 NVIDIA GeForce RTX 4090 GPUs.

### Obtaining posterior probability with the omniscient neural network model

To estimate the posterior probability of each set of rate constants, we
posited a fundamental assumption: the cells are entrapped randomly within the
microfluidic chip’s cavities. The fluorescence emissions emanating from these
cells follow a distribution dictated by the stochastic dynamics inherent to
the optical-controlled gene expression system. The trained neural network
model provides an accurate approximation for the state joint distribution
given the rate constants and the time value
$\hat{p}^{\boldsymbol{\theta}}(\boldsymbol{x}_{t}|\boldsymbol{\sigma},\boldsymbol{x}_{0},t_{0})$,
$p(\boldsymbol{x}|\boldsymbol{\sigma})$ for short. For parameter inferencing,
we used the Metropolis-Hastings algorithm with the state likelihood provided
by the neural network model. At Monte Carlo step $s$, we used a Gaussian
kernel for generating the next step proposal rate parameters
$\boldsymbol{\sigma}_{s+1}=\boldsymbol{\sigma}_{s}+\boldsymbol{l}\cdot\mathcal{N}(0,\Sigma)$,
where covariance matrix $\Sigma$ is diagonal for controlling the Monte Carlo
step size, and $\boldsymbol{l}$ is the scale factor. Since the symmetry of
Gaussian kernels, it is evident that
$p(\boldsymbol{\sigma}_{s+1}|\boldsymbol{\sigma}_{s})=p(\boldsymbol{\sigma}_{s}|\boldsymbol{\sigma}_{s+1})$.
We also use a Gamma probability density function to account for the prior
knowledge of $\boldsymbol{\sigma}$. The accepting criterion for Metropolis-
Hastings algorithm is then

$\gamma=\frac{p(\boldsymbol{\sigma}_{s+1}|\boldsymbol{x})}{p(\boldsymbol{\sigma}_{s}|\boldsymbol{x})}\cdot\frac{p(\boldsymbol{\sigma}_{s}|\boldsymbol{\sigma}_{s+1})}{p(\boldsymbol{\sigma}_{s+1}|\boldsymbol{\sigma}_{s})}=\frac{p(\boldsymbol{x}|\boldsymbol{\sigma}_{s+1})\Gamma(\boldsymbol{\sigma}_{s+1})}{p(\boldsymbol{x}|\boldsymbol{\sigma}_{s})\Gamma(\boldsymbol{\sigma}_{s})}$

which is the product of the ratio of the joint likelihoods and the ratio of
the prior probability. The movement from the original rate parameters
$\boldsymbol{\sigma}_{s}$ to the new rate parameters
$\boldsymbol{\sigma}_{s+1}$ is only allowed with the probability
$\min\left(1,\gamma\right)$. During the Monte Carlo sampling, $n$ states are
randomly chosen from the observed trajectories that correspond to different
times. The accepting criterion can then be written as:

$\gamma=\prod_{i}^{n}\frac{p(x_{i}|\boldsymbol{\sigma}_{s+1})}{p(x_{i}|\boldsymbol{\sigma}_{s})}\frac{\Gamma(\boldsymbol{\sigma}_{s+1})}{\Gamma(\boldsymbol{\sigma}_{s})}$

The trained neural network model was then used to approximate
$p(x_{i}|\boldsymbol{\sigma}_{s},t_{i})$ for estimating $\gamma$:

$\gamma\approx\exp\left[\sum_{i}^{n}\left(\ln\hat{p}_{\boldsymbol{\theta}}(x_{i}|\boldsymbol{\sigma}_{s+1},t_{i})-\ln\hat{p}_{\boldsymbol{\theta}}(x_{i}|\boldsymbol{\sigma}_{s},t_{i})\right)\right]\cdot\frac{\Gamma(\boldsymbol{\sigma}_{s+1})}{\Gamma(\boldsymbol{\sigma}_{s})}$

## Data & code availability

The met package used in this work is available in the git repositories:

https://github.com/cbliu111/cell_seg_track_software

and

https://github.com/cbliu111/met2

## Acknowledgments

C.-B.L. thanks the supports from the National Natural Science Foundation of
China Grant 32000888, the Scientific Instrument Developing Project of the
Chinese Academy of Sciences Grant YJKYYQ20180038, Jilin Province Science and
Technology Development Plan Grant 20230101152JC.

## Author contributions

E.L.E and C.B.L. designed the research and contributed to the experimental
materials. Y.F., L.L. and C.B.L. constructed the yeast strains. C.B.L.
performed the experiments, implemented the image analysis software and the
MET2 package, analyzed the data, performed the numerical experiments, and
wrote the original draft. J.W. provided critical suggestions, revised the
manuscript, and supervised the study.

## Competing interests

The authors declare no competing interests.

## Correspondence

Correspondence and requests for materials should be addressed to Jin Wang.

## SI Appendix

Supplementary information is provided in the SI.pdf file.

## References

  * 1 Elmentaite, R., Domínguez Conde, C., Yang, L. & Teichmann, S. A.  Single-cell atlases: Shared and tissue-specific cell types across human organs.  _Nature Reviews Genetics_ 23, 395–410 (2022).
  * 2 Hilfinger, A. & Paulsson, J.  Separating intrinsic from extrinsic fluctuations in dynamic biological systems.  _Proceedings of the National Academy of Sciences_ 108, 12167–12172 (2011).
  * 3 Fang, X., Kruse, K., Lu, T. & Wang, J.  Nonequilibrium physics in biology.  _Reviews of Modern Physics_ 91, 045004 (2019).
  * 4 Inoue, K. _et al._ Oscillation dynamics underlie functional switching of nf-$\kappa$b for b-cell activation.  _npj Systems Biology and Applications_ 2, 16024 (2016).
  * 5 Purvis, J. E. _et al._ P53 dynamics control cell fate.  _Science_ 336, 1440–1444 (2012).
  * 6 Gupta, A. _et al._ Inferring gene regulation from stochastic transcriptional variation across single cells at steady state.  _Proceedings of the National Academy of Sciences_ 119, e2207392119 (2022).
  * 7 Kiviet, D. J. _et al._ Stochasticity of metabolism and growth at the single-cell level.  _Nature_ 514, 376–379 (2014).
  * 8 Labib, M. & Kelley, S. O.  Single-cell analysis targeting the proteome.  _Nature Reviews Chemistry_ 4, 143–158 (2020).
  * 9 Munsky, B., Neuert, G. & Van Oudenaarden, A.  Using gene expression noise to understand gene regulation.  _Science_ 336, 183–187 (2012).
  * 10 Boersma, S. _et al._ Multi-color single-molecule imaging uncovers extensive heterogeneity in mrna decoding.  _Cell_ 178, 458–472.e19 (2019).
  * 11 Trcek, T., Lionnet, T., Shroff, H. & Lehmann, R.  mrna quantification using single-molecule fish in drosophila embryos.  _Nature Protocols_ 12, 1326–1348 (2017).
  * 12 Boyden, E. S., Zhang, F., Bamberg, E., Nagel, G. & Deisseroth, K.  Millisecond-timescale, genetically targeted optical control of neural activity.  _Nature Neuroscience_ 8, 1263–1268 (2005).
  * 13 Emiliani, V. _et al._ Optogenetics for light control of biological systems.  _Nature Reviews Methods Primers_ 2, 55 (2022).
  * 14 Nguyen, N. T. _et al._ Nano-optogenetic engineering of car t cells for precision immunotherapy with enhanced safety.  _Nature Nanotechnology_ 16, 1424–1434 (2021).
  * 15 Tan, P., He, L., Huang, Y. & Zhou, Y.  Optophysiology: Illuminating cell physiology with optogenetics.  _Physiological Reviews_ 102, 1263–1325 (2022).
  * 16 Zhao, E. M. _et al._ Optogenetic regulation of engineered cellular metabolism for microbial chemical production.  _Nature_ 555, 683–687 (2018).
  * 17 Browning, D. F. & Busby, S. J. W.  Local and global regulation of transcription initiation in bacteria.  _Nature Reviews Microbiology_ 14, 638–650 (2016).
  * 18 Sigova, A. A. _et al._ Transcription factor trapping by rna in gene regulatory elements.  _Science_ 350, 978–981 (2015).
  * 19 Schier, A. C. & Taatjes, D. J.  Structure and mechanism of the rna polymerase ii transcription machinery.  _Genes & Development_ 34, 465–488 (2020).
  * 20 Jones, D. L., Brewster, R. C. & Phillips, R.  Promoter architecture dictates cell-to-cell variability in gene expression.  _Science_ 346, 1533–1536 (2014).
  * 21 Abner, K., Aaviksaar, T., Adamberg, K. & Vilu, R.  Single-cell model of prokaryotic cell cycle.  _Journal of Theoretical Biology_ 341, 78–87 (2014).
  * 22 Kamarthapu, V. _et al._ ppgpp couples transcription to dna repair in _e. coli_.  _Science_ 352, 993–996 (2016).
  * 23 Bintu, L. _et al._ Dynamics of epigenetic regulation at the single-cell level.  _Science_ 351, 720–724 (2016).
  * 24 Orphanides, G. & Reinberg, D.  A unified theory of gene expression.  _Cell_ 108, 439–451 (2002).
  * 25 Sohmen, D. _et al._ Structure of the bacillus subtilis 70s ribosome reveals the basis for species-specific stalling.  _Nature Communications_ 6, 6941 (2015).
  * 26 Louder, R. K. _et al._ Structure of promoter-bound tfiid and model of human pre-initiation complex assembly.  _Nature_ 531, 604–609 (2016).
  * 27 Alberts, B.  _Molecular Biology of the Cell_ (CRC Press, an imprint of Garland Science, Boca Raton, FL, 2017), sixth edition edn.
  * 28 Cramer, P.  Organization and regulation of gene transcription.  _Nature_ 573, 45–54 (2019).
  * 29 Chen, J., Boyaci, H. & Campbell, E. A.  Diverse and unified mechanisms of transcription initiation in bacteria.  _Nature Reviews Microbiology_ 19, 95–109 (2021).
  * 30 Girbig, M., Misiaszek, A. D. & Müller, C. W.  Structural insights into nuclear transcription by eukaryotic dna-dependent rna polymerases.  _Nature Reviews Molecular Cell Biology_ 23, 603–622 (2022).
  * 31 Ay, A. & Arnosti, D. N.  Mathematical modeling of gene expression: A guide for the perplexed biologist.  _Critical Reviews in Biochemistry and Molecular Biology_ 46, 137–151 (2011).
  * 32 Rué, P. & Garcia-Ojalvo, J.  Modeling gene expression in time and space.  _Annual Review of Biophysics_ 42, 605–627 (2013).
  * 33 Shahrezaei, V. & Swain, P. S.  Analytical distributions for stochastic gene expression.  _Proceedings of the National Academy of Sciences_ 105, 17256–17261 (2008).
  * 34 Elowitz, M. B., Levine, A. J., Siggia, E. D. & Swain, P. S.  Stochastic gene expression in a single cell.  _Science_ 297, 1183–1186 (2002).
  * 35 Thattai, M. & van Oudenaarden, A.  Intrinsic noise in gene regulatory networks.  _Proceedings of the National Academy of Sciences_ 98, 8614–8619 (2001).
  * 36 Beaumont, M. A.  Approximate bayesian computation.  _Annual Review of Statistics and Its Application_ 6, 379–403 (2019).
  * 37 Gillespie, D. T.  Exact stochastic simulation of coupled chemical reactions.  _The Journal of Physical Chemistry_ 81, 2340–2361 (1977).
  * 38 Kilic, Z., Schweiger, M., Moyer, C., Shepherd, D. & Pressé, S.  Gene expression model inference from snapshot rna data using bayesian non-parametrics.  _Nature Computational Science_ 3, 174–183 (2023).
  * 39 Nash, A. I. _et al._ Structural basis of photosensitivity in a bacterial light-oxygen-voltage/helix-turn-helix (lov-hth) dna-binding protein.  _Proceedings of the National Academy of Sciences_ 108, 9449–9454 (2011).
  * 40 Lee, J. B., Caywood, L. M., Lo, J. Y., Levering, N. & Keung, A. J.  Mapping the dynamic transfer functions of eukaryotic gene regulation.  _Cell Systems_ 12, 1079–1093.e6 (2021).
  * 41 Shimizu-Sato, S., Huq, E., Tepperman, J. M. & Quail, P. H.  A light-switchable gene promoter system.  _Nature Biotechnology_ 20, 1041–1044 (2002).
  * 42 Li, H., Burgie, E. S., Gannam, Z. T. K., Li, H. & Vierstra, R. D.  Plant phytochrome b is an asymmetric dimer with unique signalling potential.  _Nature_ 604, 127–133 (2022).
  * 43 Smith, R. W. _et al._ Interactions between phyb and pif proteins alter thermal reversion reactions _in vitro_.  _Photochemistry and Photobiology_ 93, 1525–1531 (2017).
  * 44 Klose, C. _et al._ Systematic analysis of how phytochrome b dimerization determines its specificity.  _Nature Plants_ 1, 15090 (2015).
  * 45 Castillon, A., Shen, H. & Huq, E.  Phytochrome interacting factors: Central players in phytochrome-mediated light signaling networks.  _Trends in Plant Science_ 12, 514–521 (2007).
  * 46 Smith, H. & Jackson, G. M.  Rapid phytochrome regulation of wheat seedling extension: Light pretreatment extends coupling time, increases response lag, and decreases light sensitivity.  _Plant Physiology_ 84, 1059–1062 (1987).
  * 47 Khmelinskii, A. _et al._ Tandem fluorescent protein timers for in vivo analysis of protein dynamics.  _Nature Biotechnology_ 30, 708–714 (2012).
  * 48 Motta-Mena, L. B. _et al._ An optogenetic gene expression system with rapid activation and deactivation kinetics.  _Nature Chemical Biology_ 10, 196–202 (2014).
  * 49 Benzinger, D. & Khammash, M.  Pulsatile inputs achieve tunable attenuation of gene expression variability and graded multi-gene regulation.  _Nature Communications_ 9, 3521 (2018).
  * 50 Li, F.-W. _et al._ Phytochrome diversity in green plants and the origin of canonical plant phytochromes.  _Nature Communications_ 6, 7852 (2015).
  * 51 Legris, M., Ince, Y. Ç. & Fankhauser, C.  Molecular mechanisms underlying phytochrome-controlled morphogenesis in plants.  _Nature Communications_ 10, 5219 (2019).
  * 52 Klose, C., Nagy, F. & Schäfer, E.  Thermal reversion of plant phytochromes.  _Molecular Plant_ 13, 386–397 (2020).
  * 53 Matsushita, T., Mochizuki, N. & Nagatani, A.  Dimers of the n-terminal domain of phytochrome b are functional in the nucleus.  _Nature_ 424, 571–574 (2003).
  * 54 Burgie, E. S. _et al._ Differing biophysical properties underpin the unique signaling potentials within the plant phytochrome photoreceptor families.  _Proceedings of the National Academy of Sciences_ 118, e2105649118 (2021).
  * 55 Chen, E., Lapko, V. N., Lewis, J. W., Song, P.-S. & Kliger, D. S.  Mechanism of native oat phytochrome photoreversion: A time-resolved absorption investigation.  _Biochemistry_ 35, 843–850 (1996).
  * 56 Björling, A. _et al._ Structural photoactivation of a full-length bacterial phytochrome.  _Science Advances_ 2, e1600920 (2016).
  * 57 Sineshchekov, V. A.  Evidence for the existence of two phytochrome a populations.  _Journal of Photochemistry and Photobiology B: Biology_ 28, 53–55 (1995).
  * 58 Schmidt, P. _et al._ The complexity of the p r to p fr phototransformation kinetics is an intrinsic property of native phytochrome*.  _Photochemistry and Photobiology_ 68, 754–761 (1998).
  * 59 Burgie, E. S. _et al._ Crystallographic and electron microscopic analyses of a bacterial phytochrome reveal local and global rearrangements during photoconversion.  _Journal of Biological Chemistry_ 289, 24573–24587 (2014).
  * 60 Takala, H., Björling, A., Linna, M., Westenhoff, S. & Ihalainen, J. A.  Light-induced changes in the dimerization interface of bacteriophytochromes.  _Journal of Biological Chemistry_ 290, 16383–16392 (2015).
  * 61 Hennig, L. & Schäfer, E.  Both subunits of the dimeric plant photoreceptor phytochrome require chromophore for stability of the far-red light-absorbing form.  _Journal of Biological Chemistry_ 276, 7913–7918 (2001).
  * 62 Sweere, U. _et al._ Interaction of the response regulator arr4 with phytochrome b in modulating red light signaling.  _Science_ 294, 1108–1111 (2001).
  * 63 Ni, M., Tepperman, J. M. & Quail, P. H.  Pif3, a phytochrome-interacting factor necessary for normal photoinduced signal transduction, is a novel basic helix-loop-helix protein.  _Cell_ 95, 657–667 (1998).
  * 64 Golonka, D. _et al._ Deconstructing and repurposing the light-regulated interplay between arabidopsis phytochromes and interacting factors.  _Communications Biology_ 2, 448 (2019).
  * 65 Nagano, S. _et al._ Structural insights into photoactivation and signalling in plant phytochromes.  _Nature Plants_ 6, 581–588 (2020).
  * 66 Yoo, C. Y. _et al._ Direct photoresponsive inhibition of a p53-like transcription activation domain in pif3 by arabidopsis phytochrome b.  _Nature Communications_ 12, 5614 (2021).
  * 67 Weber, M. F. & Frey, E.  Master equations and the theory of stochastic path integrals.  _Reports on Progress in Physics_ 80, 046601 (2017).
  * 68 Cao, Y. & Liang, J.  Adaptively biased sequential importance sampling for rare events in reaction networks with comparison to exact solutions from finite buffer dcme method.  _The Journal of Chemical Physics_ 139, 025101 (2013).
  * 69 Pelechano, V., Chávez, S. & Pérez-Ortín, J. E.  A complete set of nascent transcription rates for yeast genes.  _PLoS ONE_ 5, e15442 (2010).
  * 70 Schwanhäusser, B. _et al._ Global quantification of mammalian gene expression control.  _Nature_ 473, 337–342 (2011).
  * 71 Sharma, A. K. _et al._ A chemical kinetic basis for measuring translation initiation and elongation rates from ribosome profiling data.  _PLOS Computational Biology_ 15, e1007070 (2019).
  * 72 Manavalan, B., Kuwajima, K. & Lee, J.  Pfdb: A standardized protein folding database with temperature correction.  _Scientific Reports_ 9, 1588 (2019).
  * 73 Guerra, P., Vuillemenot, L.-A., Rae, B., Ladyhina, V. & Milias-Argeitis, A.  Systematic _in vivo_ characterization of fluorescent protein maturation in budding yeast.  _ACS Synthetic Biology_ 11, 1129–1141 (2022).
  * 74 Timms, R. T. _et al._ A glycine-specific n-degron pathway mediates the quality control of protein _n_ -myristoylation.  _Science_ 365, eaaw4912 (2019).
  * 75 Varshavsky, A.  N-degron and c-degron pathways of protein degradation.  _Proceedings of the National Academy of Sciences_ 116, 358–366 (2019).
  * 76 Gligorovski, V., Sadeghi, A. & Rahi, S. J.  Multidimensional characterization of inducible promoters and a highly light-sensitive lov-transcription factor.  _Nature Communications_ 14, 3810 (2023).
  * 77 Taniguchi, Y. _et al._ Quantifying _e. coli_ proteome and transcriptome with single-molecule sensitivity in single cells.  _Science_ 329, 533–538 (2010).
  * 78 Bar-Even, A. _et al._ Noise in protein expression scales with natural protein abundance.  _Nature Genetics_ 38, 636–643 (2006).
  * 79 Shahrezaei, V. & Swain, P. S.  The stochastic nature of biochemical networks.  _Current Opinion in Biotechnology_ 19, 369–374 (2008).
  * 80 Shao, J. _et al._ Smartphone-controlled optogenetically engineered cells enable semiautomatic glucose homeostasis in diabetic mice.  _Science Translational Medicine_ 9, eaal2298 (2017).
  * 81 Pédelacq, J.-D., Cabantous, S., Tran, T., Terwilliger, T. C. & Waldo, G. S.  Engineering and characterization of a superfolder green fluorescent protein.  _Nature Biotechnology_ 24, 79–88 (2006).
  * 82 Li, Y. _et al._ A programmable fate decision landscape underlies single-cell aging in yeast.  _Science_ 369, 325–329 (2020).
  * 83 Kirillov, A. _et al._ Segment anything (2023).  2304.02643.

Figure 1: Perturbing gene expression with light. (A) Schematic of the PhyB-
based optical-controlled gene expression system. (B) Representative time-lapse
images of light-activated yeast cell incubated in traps of microfluidic chip
channels. Signal of sfGFP expression is shown with green colors. The capturing
times of the images is labeled on top of these images with unit minutes. (C)
Mean fluorescence intensity of cells as function of observation time.
Intensity of each cell is plotted as light-gray dotted lines, mean value of
these trajectories is denoted by the solid blue line and light blue region
indicates the standard deviation. Red light activation is represented by red
pulse starting at 120 min and lasts 2 mins. Far-red light deactivation is
represented by light-red box starting at 422 min and lasts 490 mins.  Figure
2: Minimized model for photoconversion and thermal reversion dynamics of
phytochromes. (A) Schematic of the state conversion model. Pr and Pfr state
are labeled with light and dark red, respectively. Fragment of PIF was denoted
as Pa. In the minimized model, dimer photoconversion and transformation
between Pfr1 and Pfr2 state are omitted. (B) Fitting of stochastic reaction
networks to the photoconversion processes. Experimental data representing the
photoconversion of PhyB(NT) at 25 ∘C (empty circle) were obtained from Figure
2 of 54. (C) Fitting of stochastic reaction networks to the thermal reversion
processes. Experimental data representing the thermal reversion of PhyB(NT)
under different concentration of PIFs at 22 ∘C (empty circle) were obtained
from Figure 3 of 43. The system was photo-activated under red light for 10 min
before darkness for observation. The photoconversion rates obtained from (B)
was used in the fitting of thermal reversion rates since the photoconversion
rates are not sensitive to temperature alterations 54, 52.  Figure 3:
Stochastic modeling and parameter inference of optical-controlled gene
expression system. (A) Stochastic modeling of the PhyB-based optical-
controlled gene expression system. Gene, mRNA, immature protein are labeled
with capital letters “G”, “R” and “U” for each species, respectively. PhyB,
PIF and fluorescence reporter are represented by letter “r”, “a” and “R”,
respectively. e.g. GR represents the gene of fluorescence reporter, and
GRPfr1Pa is the trimmer comprised by gene of fluorescence report, the
activated Pfr1 and PIF3. Pr and Pfr state are labeled with light and dark red,
respectively, as in Fig. 2. Fluorescence reporter and of PIF was denoted as
Pa. The reporter gene (GR) In the minimized model, dimer photoconversion and
transformation between Pfr1 and Pfr2 state are omitted. (B) Validation of the
trained neural network for prediction of population mean and standard
deviations of the reporter protein for 2,400 parameter combinations. Means and
standard deviations are calculated using 1,000 samples from either neural
network predictions or stochastic simulation results. (C) Posterior
probability distributions obtained through Markov Chain Monte Carlo sampling
for $\sigma_{PR}$, $\alpha_{PrPa}$, $\beta_{mRNA}$ and $\beta_{PR}$. Prior
distributions are shown by dashed dark red lines. (D) Histograms of predicted
reporter protein number compared to experimental results.  Figure 4:
Modulation of fluctuations with optical-controlled gene expression system. (A)
Steady-state probability distributions corresponding to different activation
light intensities. (B) Steady-state probability distributions of wild-type and
N-terminal mutagenesis fluorescent reporter. The inset figures in (A) and (B)
showed the change of mean fluorescence intensities along with observation
times. (C) Dependency of rate constants on activation light intensities. The
inferencing rate constants are represented as solid dots, and dash-dot lines
show the smoothed spline interpolation of the data. (D) Dependency of
logarithm of coefficient of variations (CV) on activation and degradation rate
constants (ticks are shown in lower and upper x-axes). The solid dots
represent the results obtained from stochastic simulations, and colored stars
denote the experimental results. The red solid line indicates a linear fitting
of the rate constant (Pr$\to$Pfr2) to the logarithm of CV and blue solid line
indicates an exponentially fitting of the degradation rate constant to
logarithm of CV. The red dash-dot line is an exponentially fitting of the rate
constant to the logarithm of CV. Each experimental result represents a
collection of mean fluorescent intensities from about 20,000 cells.