Noise suppression in stochastic genetic circuits using PID controllers

Inside individual cells, protein population counts are subject to molecular noise due to low copy numbers and the inherent probabilistic nature of biochemical processes. We investigate the effectiveness of proportional, integral and derivative (PID) based feedback controllers to suppress protein count fluctuations originating from two noise sources: bursty expression of the protein, and external disturbance in protein synthesis. Designs of biochemical reactions that function as PID controllers are discussed, with particular focus on individual controllers separately, and the corresponding closed-loop system is analyzed for stochastic controller realizations. Our results show that proportional controllers are effective in buffering protein copy number fluctuations from both noise sources, but this noise suppression comes at the cost of reduced static sensitivity of the output to the input signal. In contrast, integral feedback has no effect on the protein noise level from stochastic expression, but significantly minimizes the impact of external disturbances, particularly when the disturbance comes at low frequencies. Counter-intuitively, integral feedback is found to amplify external disturbances at intermediate frequencies. Next, we discuss the design of a coupled feedforward-feedback biochemical circuit that approximately functions as a derivate controller. Analysis using both analytical methods and Monte Carlo simulations reveals that this derivative controller effectively buffers output fluctuations from bursty stochastic expression, while maintaining the static input-output sensitivity of the open-loop system. In summary, this study provides a systematic stochastic analysis of biochemical controllers, and paves the way for their synthetic design and implementation to minimize deleterious fluctuations in gene product levels.


I. Comparison of analytical results with stochastic simulations
We perform stochastic simulations for the genetic circuits with different controllers presented in the main paper, using the Gillespie algorithm [1]. The analytical results for the target protein (Y ) noise obtained using the linear noise approximation are compared with the numerical simulations. For this, we compute the noise in the target protein (Y ) at the steady-state (from a large number of samples, ∼ 10 6 ) as a function of the feedback gain, keeping the mean levels constant.

Strategy for changing feedback gain maintaining constant mean levels
For the proportional and integral controller implementations, the repression from the sensor (Z) to the target protein Y is incorporated via a Hill function, where h is the Hill coefficient and z c is the amount of Z at which half-maximal repression is achieved. Then, the feedback gain given by Eq. (25) and Eq. (40) (main manuscript) reduces to dg(z) dz z= z = h(1 − g( z )).
To increase feedback gain, we vary h value. The simulation results presented below are obtained by choosing z c = z , where g( z ) = 1/2 and f p , f i = h/2. As g( z ) is independent of h, it follows from Eq. (21), (23) and (41) in the main manuscript, that the steady state mean levels are independent of the feedback gain.

Simulation results for the proportional controller
For all the simulation results presented in this SI, we consider shifted geometric distributions for the production of species: where, B is the average burst size. In Fig A,   The noise in Y copy number (CV 2 Y ) is plotted against the feedback gain for different values of the sensor noise (CV 2 Z = 0.02, 0.04, and 0.08). The lines represent analytical estimates and symbols stand for corresponding results from stochastic simulations. The simulations results show good agreement with their analytical estimate for a small sensor noise. For a high sensor noise, the quantitative agreement is poor for higher feedback gains. For this plot, the steady-state mean values of copy numbers in X, Y , and Z are kept fixed to x = 100, y = 400, and z = 100. Feedback gain is changed by varying h. Other parameter values are: degradation rates -γ x = 1.0, γ y = 3.0, γ z = 15.0; average burst sizes -B x = 2, and B y = 8. We vary B z to increase the sensor noise.

Implementation of zeroth order decay for the integral controller circuit
For this, we assume the Michaelis-Menten decay kinetics [2] for the sensor protein where v max is the maximum decay rate achieved by the enzyme-substrate interaction and k M is the amount of Z needed to attain half of the maximum decay rate. In the limit of k M z, the decay follows the zeroth-order kinetics. In this limit, the deterministic equation for the Z dynamics satisfy Eq. (36) (main manuscript) for the following choice of v max and k M : Simulation results for the integral controller In Fig B, we compare analytical and simulation results for the integral controller. The quantitative agreement between them is good for smaller values of the sensor noise. As in the case of the proportional controller, a deviation is observed for a high sensor noise with larger values of the feedback gain. The qualitative behavior, however, follows to the analytical prediction. This deviation is attributed to the linearization approximation on the basis small noise assumption which violates in the presence of large fluctuations and strong nonlinearity. The noise in Y copy number (CV 2 Y ) is plotted against the feedback gain for different values of the sensor noise (CV 2 Z = 0.01, 0.02, and 0.03). The lines represent analytical estimates and symbols stand for corresponding results from stochastic simulations. The simulations results show good agreement with their analytical estimate for a small sensor noise. For a high sensor noise, the quantitative agreement is poor for higher feedback gains. For this plot, the steady-state mean values of copy numbers in X, Y , and Z are kept fixed to x = 100, y = 200, and z = 100. Feedback gain is changed by varying h. Other parameter values are: γ x = 1.0, γ y = 3.0, k z = γ y , B x = B y = 2. We vary B z to increase the sensor noise.

II. A combined circuit of proportional and integral controllers
In the main paper, we consider the proportional, integral, and derivative controller separately. How does a combined PI (proportional and integral) implementation work in terms of target protein noise as opposed to the individual controllers? To study this, we consider a combined PI system in the presence of extrinsic noise and solve numerically. The circuit is schematically shown in Fig C(i). This system consists of four species X (for extrinsic noise), Y (target protein), Z p (proportional controller), and Z i (integral controller). The species Y activates both controllers linearly. Both controller species Z p and Z i repress the synthesis of the target protein Y , and this repression is modelled the following two-dimensional function, where h p and h i are the Hill coefficients and z c is the value of repression coefficient which we assume to be the same for both the species for simplicity [3]. The dynamics of external noise is implemented in the same way as described in the main text. The stochastic dynamics of Z p , Z i , and Y are given by, Note that the zero-order decay of Z i is implemented using Michaelis-Menten decay kinetics as described above (Eqs. S4 and S5). For simplicity, we choose the burst size distributions for both controllers identical and independent. For the combined repression function (Eq. S6), the proportional and integral feedback gain can be defined as, To increase feedback gains, we vary the value of Hill coefficients h p and h i . The simulation results presented below are obtained by choosing We simulate the stochastic dynamics of X, Y , Z p , and Z i using Gillespie algorithm and compute the noise in the target proteins at the steady-state from a large number of sample points. In Fig C, we plot the target protein noise (normalized by target protein noise for open loop circuit) as a function of the proportional feedback gain f p and integral feedback gain f i . When the extrinsic noise is small, the only proportional controller case (f i = 0 in Fig C) is the best for noise reduction in Y . Interestingly, when extrinsic noise is large, the combined PI controller performs better as opposed to individual controllers (f i = 0 or f p = 0 axes).

III. Parametric constraints for derivative control
The burst frequency function for target protein synthesis, which depends on both the target protein (Y ) and controller species (Z), must obey particular constraints to act as a derivative controller. Here we discuss this limitation. Let g(y, z) be the general burst frequency function. We can linearize it around the mean as g(y, z) − g y , z ≈ g y , z S z z z + S y y y .
(S10)  Here, S z and S y are normalized sensitivities of g(y, z) with respect to Z and Y respectively given by (S11) Using Eq. (49) in the main text, and assuming γ z is large (i.e., fast sensor dynamics), the Laplace transform of the right-hand-side of (S10) is Here, the first term on the RHS represents proportional control while the second term is corresponds to derivative control. Therefore, this represents combined proportional derivative controller. A pure derivative controller is obtained only when S z = −S y for which the first term drops out of the equation.
For the regulation of burst frequency Y as shown in the main text, with Hill functions for the activation and repression we get Eq. (50) in the main text. Using S z = −S y derived above and (S11) we get the following relationship This relationship simplifies to h z = h y = h for z c z and y c y and leads to g(y, z) = (z/y) h as reported in the main text.
In this analysis we show that while the general derivative control requires S z = −S y , for simplicity we choose the formulation consisting of Hill functions with strong binding affinity of the regulators with the same Hill coefficients.

IV. Dependence of target protein noise on Y → Z activation strength
In the main text, we mostly discussed the results for a given activation strength of Y → Z reaction, i.e. for a fixed z value. Here, we analyze how the target protein noise depends on the z .
In the proportional controller, we can change z by varying k z which represents the production rate of Z due to the activation by Y given by Eq. (21) in the main manuscript. This changes CV 2 Z according to Eq. (31) and changes f p based on Eq. (27). For a fair comparison of CV 2 Y with changing k z we need to keep y constant. To keep y constant, we vary the maximum transcription rate k y following Eq. (23). With increasing z , the feedback strength increases according to So, for low z , f p ∝ h z h and given that the sensor noise component for low z becomes Sensor noise  Fig D (i,ii). Note that for higher h, this effect is amplified due to the presence of the scaling of the sensor noise component with respect to h. The other two components are monotonically buffered with increase in z which reflects the sole increase in f p in For the integral controller, z is a parameter which can be increased independent of k z due to the way we have defined its dynamics given in Eq. (36) in the main manuscript. CV 2 Z and f i change in a similar fashion to those in proportional controller. We compensate with k y to keep y constant based on Eq. (41). The sensor noise in this case given by This shows that for h = 1, the noise decreases, while for h > 1, it follows the same profile as in the proportional controller case as seen in Fig D (iii,iv). While the noise from stochastic expression is independent of f i and CV 2 Z , the extrinsic noise follows which is an increasing function with respect to z for the regime where h = 1 (Fig D (iii)). In the highly non-linear regime h = 6, f i increases steeply at intermediate values of z and causes a slight U-shaped profile as seen in Fig D (iv). For the derivative controller, we can change z by varying k z which represents the production rate of Z due to the activation by Y given by Eq. (48) in the main manuscript. To keep y constant, we vary the maximum transcription rate k y following Eq. (53). In this analysis, the feedback gain f d is independent of z and we find that CV 2 Y follows the change in CV 2 Z . CV 2 Z reduces asymptotically with increasing level of Z as seen in (S16). This gives the monotonically decreasing profile as seen in Fig D (v,vi).

V. External disturbance as an OU process
In this section we derive the expression for CV 2 Y without feedback using the dynamics of X given by: Here, x is the mean level of the external disturbance, γ x is the time scale of the disturbance, σ x represents strength of the disturbance and w(t) is the Wiener process. The hybrid system is then given by the above OU process along with the dynamics of Y shown in Eq. (13) in the main text. The resulting moment dynamics for an arbitrary function φ(x, y) using an extended generator L of   The noise in Y for the derivative controller with respect to z for h = 1 and h = 6 respectively is shown. Parameters used are: CV 2 int = 0.25, CV 2 X = 0.7, γ z = 3γ y = 5γ x , and B z = 2 with shifted geometrically distributed bursts. Noise in Y is obtained from Eq. (57).
The resulting second order moment equations are found by appropriately substituting for the function φ(x, y) y 2 dt = −2γ y y 2 + γ y y + k y B 2 y + 2 B y k y xy x (S24) xy dt = γ x ( x y − xy ) + k y B y x 2 x − γ y xy (S25) Solving these equations and obtaining the expression for CV 2 Y in terms of CV 2 X we see that is identical to Eq. (17) in the main text where CV 2 X = σ 2 x 2γx x 2 instead of Eq. (12) in the main text.