Effects of spatial heterogeneity on bacterial genetic circuits

Intracellular spatial heterogeneity is frequently observed in bacteria, where the chromosome occupies part of the cell’s volume and a circuit’s DNA often localizes within the cell. How this heterogeneity affects core processes and genetic circuits is still poorly understood. In fact, commonly used ordinary differential equation (ODE) models of genetic circuits assume a well-mixed ensemble of molecules and, as such, do not capture spatial aspects. Reaction-diffusion partial differential equation (PDE) models have been only occasionally used since they are difficult to integrate and do not provide mechanistic understanding of the effects of spatial heterogeneity. In this paper, we derive a reduced ODE model that captures spatial effects, yet has the same dimension as commonly used well-mixed models. In particular, the only difference with respect to a well-mixed ODE model is that the association rate constant of binding reactions is multiplied by a coefficient, which we refer to as the binding correction factor (BCF). The BCF depends on the size of interacting molecules and on their location when fixed in space and it is equal to unity in a well-mixed ODE model. The BCF can be used to investigate how spatial heterogeneity affects the behavior of core processes and genetic circuits. Specifically, our reduced model indicates that transcription and its regulation are more effective for genes located at the cell poles than for genes located on the chromosome. The extent of these effects depends on the value of the BCF, which we found to be close to unity. For translation, the value of the BCF is always greater than unity, it increases with mRNA size, and, with biologically relevant parameters, is substantially larger than unity. Our model has broad validity, has the same dimension as a well-mixed model, yet it incorporates spatial heterogeneity. This simple-to-use model can be used to both analyze and design genetic circuits while accounting for spatial intracellular effects.

with domain D(L v ) = {y ∈ L 2 (0, 1) : v 2 d dx y v x=0,1 = 0}: Next, we introduce the Hilbert space L 2 v (0, 1) that we use for our analysis. This space is isomorphic to L 2 (0, 1) however, the operator (1) is self-adjoint with respect to the inner product in L 2 v (0, 1). Definition 2. (Weighted L 2 (0, 1) space) For smooth v : [0, 1] → R + we denote L 2 v (0, 1) as a weighted space of the square integrable functions such that f ∈ L 2 v (0, 1) if and only if : ∈ L 2 (0, 1). The inner-product is defined as f, g v := 1 0 dx, for f, g ∈ L 2 v (0, 1) andv(x) = v(x) v . Furthermore, let v = [v 1 (x), . . . , v n (x)] T , where v i (x) : [0, 1] → R + is smooth, then z ∈ L 2 v ((0, 1), R n ) if z i ∈ L 2 vi (0, 1) for i = 1, . . . , n, and the inner product in this space is defined as f , for any y ∈ L 2 v (0, 1). Thus, when performing convergence analysis we may use the norms defined in L 2 v (0, 1) and L 2 (0, 1) interchangeably. III λ 1 = 0 and ψ 1 =v(x) = v(x) v Proof: The proof of (I) and (II) follow from Sturm-Liouville theory [1]. To prove (III), we take the weighted inner product of both side of L v ψ i = λ i ψ i with ψ i , we use orthonormality, use integration by parts, and apply the boundary conditions: The maximum of (2) is achieved (λ 1 = 0) for ψ 1 =v(x), since substitutingv(x) directly into L v , one observes that L v (v(x)) = 0 and we have that ||ψ 1 || L 2 v = ψ 1 , ψ 1 v = 1, therefore λ i < 0, ∀i > 1. The following will introduce the notion of contracting dynamical systems. A dynamical system is said to be contracting within an open and connected subspace of the state space, if all trajectories starting within this region converge exponentially to each other. We provide sufficient conditions to guarantee that a dynamical system is contracting, and finally show that contracting systems have a particular robustness property. The robustness property will be exploited several times in our analysis to perform our model reduction.
Theorem 1. Letż = f (t, z) be a dynamical system in the Hilbert space H = R n where f is a smooth nonlinear function. A dynamical system is said to be contracting within an open and connected subspace of the state space χ ⊆ H, if all trajectories starting within this region converge exponentially to each other. A sufficient condition for a system to be contracting in H is the existence a uniformly positive definite matrix P (t, z) ∈ R n×n and constant λ > 0 such that where ξ is the contraction rate of the system andṖ is the total time derivative of P .
Lemma 2. (Hierarchies of contracting systems) Let be dynamical systems in the Hilbert space H = R n+m where f 1 : [0, ∞) × R n → R n and f 2 : [0, ∞) × R n × R m → R m are smooth nonlinear functions. Then, sufficient conditions for (4) to be contracting in χ = χ 1 ⊕ χ 2 , where χ 1 and χ 2 are open connected subspaces of R n and R m , respectively, are I ∂f 1 ∂z1 satisfies (3) in χ 1 for some P 1 ∈ R n×n such that P 1 = P T 1 > 0 II ∂f 2 ∂z2 satisfies (3) in χ 2 for some P 2 ∈ R m×m such that P 2 = P T 2 > 0 III ∂f 2 ∂z1 is uniformly bounded for all t ≥ 0, z 1 ∈ χ 1 , and z 2 ∈ χ 2 Proof: See hierarchal structures in [2] and applied in [3] Proof: This follows from the result of Seciton 3.7 Result vii in [2]. Let R z (0) be the length of the straight path connecting z p (0) and z(0), that is, R z (0) = zp(0) z(0) ||δz(0)|| H = ||z p (0) − z(0)|| H , where δz is the virtual displacement (see [2] for details. Each point in the straight path connecting z p (0) and z(0) evolves in time and we denote the length of the path connecting these points be given by R z (t). Precisely, R z (t) = zp z ||δz(t)|| H . Notice that ||z p (t) − z(t)|| H ≤ R z (t), (5) since the straight line path is the shortest path between the trajectories (they would be equal if the initial straight line segment remained a line for all time, but this is only the case for constant vector fields). Let R(t) = zp z ||P 1/2 δz|| H . Notice that From Equation 15 in [2], we have thatṘ + ξR ≤ ||P 1/2 d|| H and thus From (6), we have that Finally, by (5), this implies that The desired result holds for * = ζ/(2ξ) and L * = 2L 1 /ζ + L 2 /ξ.
1.2 Solutions of diffusing states converge to the null space of the differential operator L v Let v : [0, 1] → R n d + be a smooth vector-valued function and V (x) = diag(v(x)). For the state vectors z s (t, x) ∈ H s := L 2 ((0, 1), R ns ), ∀t ≥ 0 and z d (t, x) ∈ H d := L 2 v ((0, 1), R n d ), ∀t ≥ 0, consider the following reaction-diffusion system: is an in Definition 1. The elements of z d may be thought of as "freely diffusing" and those of z s as "spatially fixed". For the general system (7), we show in this section that the spatial profiles of the diffusing species z d , approach those of v(x), that is, the null space of L v when 1.
II There exists a positive constant M > 0 such that x, z s , z d ), t > 0, x ∈ (0, 1), The following theorem is a direct consequence of the operator L v from Definition 1 being self-adjoint and negative semi-definite. This result will show that the infinite dimensional left over dynamics z ⊥ d (t, x) become order after a fast transient.
Theorem 2. Consider the system defined by (7). Suppose that Assumption 1 holds and The following proof uses a similar logic as the the min-max theorem for matrices [4] to derive and upper bound for Let λ i,j and ψ i,j denote the j-th eigenvalue and eigenfunction, respectively, of . From the orthonormality of the eigenfunctions, linearity of L i v = L v i , and the ordering of eigenvalues (Lemma 1-I), notice that: Boundary conditions Here v E,i (x), v S (x), and v c,i (x), are the available volume profiles of E i , S, and c i , respectively, and L v is as in Definitions 1. The parameters D Ei , D ci , and D s , are the enzyme, complex, and substrate diffusion coefficients, respectively, is a dimensionless parameter that captures the speed of diffusion (with respect to dilution). A species being spatially fixed translates to the flux being zero throughout the whole spatial domain. In Case II, we denote the location of the fixed species E i , as x * i ∈ (0, 1) for i = 1, . . . , n. In Case III, we denote the location of the fixed species S, as x * s ∈ (0, 1).
We denoteĒ i (t),S(t), andc i (t) to be the space averaged enzyme, substrate, and complex concentrations, respectively (e.g.,Ē i (t) = 1 0 E i (t, x)dx). The dynamics governing these space averaged variables are derived by integrating (8) in space and applying the boundary conditions and are given by: where overbars denote spatially averaged variables. We make the follow assumptions necessary to state our main result Assumption 2. Consider (8), for i = 1, . . . , n, assume that I the functions α i (t, x) and α s (t, x) are smooth in each argument Table 1) are smooth, strictly greater than zero, and bounded above by unity.
The following assumption makes precise what it means for the spatially fixed species E i (in Case II) and S (in Case III) to be spatially localized at x * i and x * s , respectively. Assumption 3. (Localization of spatially fixed species) Consider the system (8). Let x * i ∈ (0, 1) for Case II and x * s ∈ (0, 1) for Case III, be given by Table 1. Let δ * = min(x 1 , . . . , x n , 1 − x 1 , . . . , 1 − x n ) for Case II and δ * = min(x * s , 1 − x * s ) for Case III. We assume that for a given δ > 0 such that δ < δ * , the functions α i (t, x) and α s (t, x) satisfy The following definition will provide the candidate reduced model to approximates (9).
Remark 3. The set Ω z , guaranteed to exist in Theorem 3, depends on δ in Cases II-III and this dependence is made precise in the proof.
Road map of proof : The rest of this section is dedicated towards proving Theorem 3. We first apply Theorem 2 to show that the spatial profile of a freely diffusing species converges to its available volume profile (e.g,. v E,i (x), v S (x), and v ci (x)) exponentially fast in the time scale associated with diffusion. If the localization assumption holds, then (9) has the form described by Definition 3, but with additional "disturbance" terms of order and δ. We proceed to demonstrate the system described in Definition 3 is contracting and apply the robustness property of contracting systems (Lemma 3) to show closeness between its solutions and those of (9).
The following result will define a positively invariant and bounded subset of R 2n+1 such that solutions to (8) starting within this set at t = 0, remain within this set for all times and spatial values. To apply Theorem 2 to (8), the existence of such a positively invariant set is required by Assumption 1.
Proof: We apply Theorem 1 in [5], which states that for a parabolic PDE system with sufficiently smooth coefficients, a closed convex subset of euclidean space is positively invariant if the vector field corresponding to the "reaction dynamics" never points outwards at the boundaries of the set. To apply this theorem we first make a coordinate transformation. The spatial differential operator (as in Definition 1) for a general diffusing species y(t, x) with available volume profile v(x), given by is not in the standard form stipulated by Equation 1.2 in [5]. Therefore, for Cases I-III, the following coordinate transformation is made With the transformation u(t, x) = y(t, x)/v(x), the differential operatorL v given bŷ is in the postulated form. Let . . , f En , f c1 , . . . , f cn , f s ] T is the vector field of the reaction dynamics of (8). The vector field corresponding to the reaction dynamics of the transformed system is given byf (t, x, u) := Λ −1 f (t, x, Λu). Let and χ u (x) = {u ∈ R 2n+1 : u i ≥ 0, ∀i and u T n u (x) ≤ Γ}. We now check that the vector fieldf (t, x, u) does not point outward at all the boundary points of χ u (x), ∀x ∈ [0, 1]. Checking thatf The set ∂χ u (x) = {u ∈ R 2n+1 : u T (x) n u (x) = Γ}, corresponds to the boundary points defined by planar surface u T (x) n u (x) = Γ with normal vector n u (x), we need to check that for all boundary point u * ∈ ∂χ u (x) we have thatf (u j, * Λ j,j + 2u j+n, * Λ j+n,j+n ) + u 2n+1, * Λ 2n+1,2n+1 ) ≤ n j=1α j +α s − ( n j=1 (u j, * + 2u j+n, * ) + u 2n+1, * )v * , Assumption 2 Thus, χ u (x) is a positively invariant set in the u(t, x) coordinates. The corresponding invariant set in the z(t, x) coordinates is given by χ.
Definition 4. Consider the systems given by (8) and (9) and let Cases I-III correspond to those in Table 1. We define For Case I and Case III: For Case I and Case II: and For Case I: Lemma 4. Consider the systems given by (8) and let > 0 be defined for Cases I-III by Table 1. Let w ⊥ (t, x) be as in Definition 4 and suppose Assumption 2 holds. Let χ ⊂ R 2n+1 be as described in Claim 1 and let z(t, x) be given by (11). Then there exists ζ, L ⊥ > 0 such that for all z(0, x) ∈ χ, ∀x ∈ [0, 1] and > 0, where m = 2n + 1 in Case I, m = 1 in Case II, and m = n in Case III.
Proof: This results follows directly from Theorem 2 where Assumption 1-I is satisfied by χ and Assumption 1-II is satisfied by the smoothness of the reaction dynamics in (8), the compactness of the sets χ and [0, 1], and the temporal boundedness (Assumption 2).
Remark 4. From the proof of Theorem 2, one can observe that L ⊥ depends on the size of χ. Once Assumption 3 is made, the size of χ will depend on δ for Case II-III. Thus, L ⊥ depends on δ for Cases II-III.
Next we define the space averaged total enzyme and substrate quantities and show that these are the same for (9) and (10). Furthermore, we will show the dynamics for these quantities are governed by uncoupled, linear, and contracting ODEs.
Remark 6. From the linear and uncoupled structure of (16), it is clear that the dynamics forĒ i,T (t) are contracting with contraction rate λ i = γ i +1 for all i = 1, . . . , n. Similarly, theS T (t) dynamics are contracting with contraction rate λ s = γ s + 1.
Remark 7. By Assumption 3, we have thatĒ * i,T andS * T are independent of δ and thus these upper bounds forĒ i (t) andS(t), respectively, are also independent δ.
The following results demonstrates that Assumption 3 implies that the concentration of fixed species is localized at the region specified in Table 1. The parameter δ > 0, controls the amount of localization. Proposition 1. Consider the systems given by (8). Let x * i ∈ (0, 1) for Case II and x * s ∈ (0, 1) for Case III, be given by Table 1. Suppose that Assumption 3 holds for x * i , x * s , and a given δ > 0. Let for Case III, we have that The following claim will aid us in rewriting (9) in the form of the reduced dynamics given by (10) with additional "disturbance" terms of order and δ. The claim is written to handle Case I-III.
1. There exists constant k 1 , k 2 > 0 such that Proof: To proof the first claim, notice that ) and leveraging the Cauchy-Schwarz inequality in H, For the second part of the claim, notice that , the existence of c is guaranteed by the mean-value theorem for integrals [6]. Notice that 1 δ, the following uses the smoothnessv 1 (x) to guarantee uniform continuity i.e., the existence of L v > 0 such that for all Finally, In the proof of Claim 3-2, y * 2 may depend on δ since we assume thatȳ * 2 is independent of δ (one expects that y * 2 increases with decreasing δ), thus k 3 may depend on δ. Corollary 2. Consider the systems given by (8). The assumptions necessary to apply Claim 3 to are satisfied by Claim 1, Proposition 1, and Claim 2 (along with Remark 7). Furthermore, considering the results from Lemma 4, we are guaranteed he existence of L i > 0 for i = 1, . . . , 5 such that for all > 0 where θ * i is given by (10b) for Cases I-III (as in Table 1), and where ζ is as in Lemma 4. The coefficients L 3 and L 4 may depend on δ by the discussion in Remark 4 and Remark 8.
is given by (10), thec(t) andĉ(t) dynamics may be written as where ∆ i = ∆ i (t) as in Corollary 2 and f c : . By the form of (22), it is clear that thec(t) dynamics (22a) are in the form theĉ(t) (22b) with additional "perturbation" terms of order and δ.
The variablesĒ T (t),S T (t), andc(t) are enough to fully describe (9) andĒ T (t),S T (t), andĉ(t) are enough to fully describe (10). Notice that and by Claim 2 these terms are uniformly bounded in time and for i = 1, . . . , n. Considering Lemma 2 with z 1 (t) = [Ē T T (t),S T (t)] T and z 2 (t) =ĉ(t), then condition (I) is satisfied by the discussion in Remark 6 and condition (III) is satisfied by (23), thus we can treatĒ T andS T (t) as time varying inputs to (22b) when showing that (II) is satisfied.
We now show that the dynamics (22b) are contracting and thus we can apply the robustness property of contracting systems (Lemma 3) to show that the solutions of (22b) and (22a) are close.
Lemma 5. Consider the system (8) and let > 0 be defined for Case I-III by Table 1. Letĉ(t) be given by (22b) andc(t) be given by (22a). Suppose the conditions of Claim 1, Lemma 4, and Claim 2 hold. Then, there exists L c,1 , * > 0, such that for all ≤ * , the solutionsĉ(t) andc(t) satisfy where |∆ c | = L c,1 for Case I. For Cases II-III, if in addition, the conditions of Proposition 1 hold for all Proof: Consider the metric , . . . , 1 a n θ * nÊn (t) whereÊ i (t) is given by (10). The total time derivative of P is given bẏ .
The Jacobian of (22b) is given by and σ is a rank one matrix given by With the chosen metric P (t), we have that since the symmetric rank one matrix σ 1 , has n − 1 zero eigenvalues and the nontrivial eigenvalue is λ σ = −n for eigenvector v n = [1, . . . , 1] T . Recalling the positivity ofÊ i (t),ĉ i (t),Ŝ(t) (Claim 1), andᾱ i (t), we have that By Theorem 1, the system (22b) is contracting with contraction rate . Therefore, for a given δ > 0, we apply the result from Lemma 3 to the nominal system (22b) and the perturbed system (22a). Let * = ζ/(2ξ) and recalling thatc i (0) =ĉ i (0), then by Lemma 3, there exists l 1 , l 2 (δ), l 3 > 0 such that for all < * where M is a constant upper bound on the square root of the ratio of the biggest and smallest eigenvalues of P (t). Thus L c,1 = M l 1 , L c,2 (δ) = M l 2 (δ), and L c,3 = M l 3 . We now show that M exists. Let r(t) be the ratio of the biggest and smallest eigenvalues of P (t). By Claim 2, there existsĒ i, * ,Ē * i,T ∈ R + independent of and δ such thatĒ i, * ≤Ê i (t) ≤Ē * i,T , ∀t ≥ 0, and thus Corollary 3. Recall Definition 5 and Remark 5, we have that where |∆ c | is as in Lemma 5. Thus, the quantity |∆ z | as claimed to exist in Theorem 3, may be given by Remark 9. We now comment on the sets Ω z ∈ R n+1 and Ωz ∈ R n+1 as claimed to exist in Theorem 3. Let χ be as in Claim 1.
In Case II, for Proposition 1 to hold we also require that for i = 1 . . . , n, that thus Ω z is the intersection between the set that satisfies this condition and χ. In Case III, for Proposition 1 to hold we also require that thus Ω z is the intersection between the set that satisfies this condition and χ. As discussed in Remark 4, χ may depend on δ.

Fast Diffusion and Bingding Dynamics
The approximation result of Theorem 3, holds well if w ⊥ is small, where w ⊥ is given by Definition 4. This was guaranteed by Lemma 4 which is based on Theorem 2. The proof of Theorem 2 was based on the principle that However, for (8), if and η i are of similar order of magnitude, and cannot be neglected. The terms and η i being comparable corresponds to diffusion and the binding between E i and S occurring at similar timescales, which often time occurs within the cell [7]. Here we show that when both diffusion and the binding dynamics dominate in (8), z ⊥ d = 0 is still the quasi-steady state, that is, all freely diffusing species mirror their available volume profile.
When both diffusion and the binding dynamics dominate in (8) we have that We compute the quasi-steady state of the "fast dynamics in (26) for Cases I-III, that is E i (t, x), c i (t, x), and S(t, x) such that Case I: In this case, (27) is satisfied for with the additional constraint that which in general (28) is a stringent condition since it is required to hold for all x ∈ [0, 1]. However, by the Case II and III: For this cases, (27) is satisfied for Thus, when both diffusion and the binding dynamics dominate, the quasi-steady states are still those that correspond to freely diffusing species converging to their available volume profile. We observed that for Case I, this was possible by the fact that v ci (x) = v Ei (x)v S (x). In a future study it should be shown that this quasi-steady state is a stable solution of the fast dynamics in (26).
2ρ . For a species diffusing inside the cell with radius of gyration r, we model the available volume profile v(x) as: v where κ is an empirical coefficient (as discussed in [8]) and r * 2 = 1/(2κπρ). From the parameter values in [8], r * ≈ 23 nm. As commented in SI Section 7, r * can be estimated for a given context (e.g., growth conditions and strain ) by analyzing the concentration profile inside the cell of a freely diffusing species with a known radius of gyration, which is possible via superresolution imaging [9]. As shown in Figure A, we estimate the chromosome density as a step function.
where v 0 = e − (r/r * ) 2 2(1−∆x) and ∆x is the distance between the end of the chromosome and the cell poles (see Figure 1 in the Main Text). Its clear now, that our choice to defineρ(x) = ρ(x) 2ρ was motivated by the fact that when ∆ x = 1/2 (nucleoid evenly spread out between mid-cell and the halfway point between mid-cell and the cell poles), we have the convenient expressionsρ(0) = 1 and v 0 = e −(r/r * ) 2 . Notice that v 0 → 0 as Bounds on θ * i : For the species E i and S as described in the Main Text with radius of gyration r e,i and r s , respectively. The available volume profiles are given by and Case I-III (as in the Main Text, Equation 15) in the Main Text, the bounds on θ * i are given by Figure A: Idealization of the chromosome density which yields a simple estimate of the available volume profile. The chromosome density is approximated as a step function implying by (30) that the available volume profile is also a step function. Here ∆x is the distance between the end of the chromosome and the cell poles.
• For Case I, this idealization implies that The upper limit of θ * i is 1 ∆x and is reached as v 0,Ei and v 0,S approach zero, which occurs as r e,i /r * → ∞ and r s /r * → ∞. The lower limit of θ * i is unity and is achieved if v 0,Ei or v 0,S approach one, which occurs if any of the two species is sufficiently small (r e,i /r * 1 or r s /r * 1 ). Since θ * i ≥ 1, it implies that the binding betwen E i and S is always equal to or greater than that predicted by a well-mixed model (Main Text Equation 9).
• For Case II-III, let x * i and x * s as in Assumption 3 in the Main Text and thus When x * i ≤ 1 − ∆x for Case II and x * s ≤ 1 − ∆x for Case III, we have that θ * i ≤ 1, the lower limit θ * i = 0 is achieved when r e /r * → ∞ for Case II (r s /r * → ∞ for Case III), the upper limit θ * i = 1 is achieved when r e /r * → 0 for Case II (r s /r * → 0 for Case III). When x * i > 1 − ∆x for Case II and x * s > 1 − ∆x for Case III, we have that θ * i ≥ 1, the lower limit θ * i = 1 is achieved when r e /r * → 0 for Case II (r s /r * → ∞ for Case III), the upper limit θ * i = 1/∆x is achieved when r e /r * → ∞ for Case II (r s /r * → ∞ for Case III).

Protein Production: Transcription and Translation
We consider gene (D) being transcribed by RNAP (S) to form a DNA-RNAP complex (c s ) to produce mRNA (m) which is translated by ribosomes (R) to form mRNA-ribosome complex (c m ) which produces protein P . The mRNA's degrade at rate γ. The RNAP, and ribosomes are produced at rates α s (t, x) and α r (t, x), respectively. We assume all species dilute at rate µ, the cells growth rate. where γ is the mRNA degradation rate, a s and d s are the association and dissociation rate constants, respectively, between RNAP and the gene D, κ s is the catalytic rate of formation of mRNA m, a m and d m are the association and dissociation rate constants, respectively, between ribosomes and mRNA , κ m is the catalytic rate of formation of protein P. We assume that the total concentration of D is conserved, so that D T (x) = D(t, x) + c s (t, x) and that D T (x) is localized at x = x * . Spatial-temporal Dynamics: The dynamics corresponding to these biochemical reactions are given by: where the spatial variable has been normalized by L (cell-length) and the time variable has been normalized by 1/µ the time scale associated with dilution. The flux dynamics and boundary conditions are given by, where v m (x), v c (x), v R (x), and v s (x) are the available volume profiles for the mRNA, mRNA-ribosome complex, ribosome, and RNAP, respectively and χ m = D m /(L 2 µ), χ c = D c /(L 2 µ), χ R = D R /(L 2 µ), and χ s = D s /(L 2 µ), are the dimensionless diffusion coefficients for the mRNA, mRNA-ribosome complex, ribosome, and RNAP, respectively. The space averaged protein concentrationP (t) is given by

Values for dimensionless parameters:
We set all production rates with respect to that of RNAP such thatᾱ S = 1 0 α S (x)dx = 1. All time scales relative to µ = 0.5 1/hr, consistent with the experiments [10] The total number of RNAP (N RNAP ) ranges between 2,000 -10,000 we took it to be 5,000 [9]. The total number of ribosomes (N ribo ) was taken to be 10, 000 and since both RNAP and ribosomes and RNAP are stable, it impliesᾱ r = 1 0 α r (x)dx = N ribo NRNAP = 2. mRNA degradation is about 10 times faster than dilution [7], therefore, γ = 10.The rate of transcription (translation) is about 80 (40) times faster than dilution [7], thus we choose κ s = 80 and κ m = 40. We assumed that the DNA is on a high copy plasmid (≈ 500 copies) and thusD T = The length of the cell is about 3µm and thus L = 1.5µm [8]. The diffusion coefficient of RNAP is taken to be D s = 0.22µm 2 /s [11] and thus χ s = 704. The diffusion coefficient of free ribosomes is taken to be D r = 0.4µm 2 /s [8] and thus χ r = 1280. In [9], the diffusion coefficient of polysomes is 0.05 ± 0.02µm 2 /s, and thus we take the diffusion coefficient of a free mRNA to be the upper bound 0.07µm 2 /s and thus χ m = χ c = 224.
For the following, the spatial profiles for the production rates are given proportional to their functional form since the constant that fully specifies them is such that the production rate per-cell satisfies the above values. Additional simulation details for Figure 5-A in the Main Text: D T (x) ∝ e −20x when DNA near mid cell and D T (x) ∝ e 20(x−1) when DNA at cell poles. The RNAP production was kept roughly spatially constant α s (x) ∝ e −.001x . The binding and unbinding coefficients for DNA-RNAP were a s = 1000 and d s = 1000. We set a m = d m = 0 such that mRNA did not bind to ribosomes and thus the free amount of mRNA is equivalent to the total mRNA. Additional simulation details for Figure 5-B in the Main Text: D T (x) ∝ e −.001x is chosen to be roughly constant. The RNAP production was kept roughly spatially constant α s (x) ∝ e −.001x . The ribosome production was kept roughly spatially constant α r (x) ∝ e −.001x . The RNAP radius of gyration was taken to be r s /r * = 0.001 such that its excluded volume effects were negligible. The binding and unbinding coefficients for DNA-RNAP were a s = 1000 and d s = 1000. The binding and unbinding coefficients for ribosome-mRNA were a m = 10 and d m = 10.  For both values of η TX the error is high at t = 0 since the initial RNAP spatial profile is chosen to be a constant but quickly decays to less than 2%. The rest of the simulation set up and parameters are identical to those of Figure 5-A.  Figure 5-A in the Main Text. When the DNA is pole localized the relative error is less than 1% and when the DNA is localized near mid-cell the error is less than 10%. where the first and second reaction model the loading and translation steps, respectively. The translation dynamics corresponding to (36) are given by   For both values of η TL the error is high at t = 0 since the initial mRNA spatial profile is chosen to be a constant but quickly decays to less than 1%. The rest of the simulation set up and parameters are identical to those of Figure 5

Multiple Ribosomes on a Single Strand of mRNA
where c s (t, x) is given by (34), L m (m) and L R (R) are given by (35) and . Integrating (37) in space yields: . (39) The production rate ofP denoted byᾱ P is given byᾱ P = κ tct (t). From our analysis in Section 1.2, we expect that R(t, x) ≈R(t)v r (x), m(t, x) ≈m(t)v m (x), and c l (t, x) ≈c l (t)v l (x) (this is verified computationally in Figure I), and thus we can estimate θ * l (t) and θ * t (t) by the constants Let K d = (d l /a l ) 1/(Nr−1) , K t = (d t +κ t )/a t , β l = (γ +1)/d l , and β t = (γ +1)/(κ t +d t ), if β l , β t , β t R/K t 1 (dilution and mRNA degradation is much slower the rate of ribosome unbinding and K t is sufficiently large), then a simple expression for the steady state protein is given bȳ wherem T =m +c l +c t = κ scs /(γ + 1) is the total mRNA.

Transcription Factor Regulation
Intracellular signaling to control gene expression is often done via transcription factors (TFs). In this section we model a general transcription factor architecture where the repressor P r dimerizes to form c 1 (e.g., TetR dimerizes before targeting gene [12]) and then blocks the transcription of gene D that produces protein P. The biochemical reactions corresponding to this process are: where α is the production rate of P r , a 1 (d 1 ) is the association (dissociation) constant to form the c 1 complex, a 2 (d 2 ) is the association (dissociation) constant to form the the c 2 complex, and κ is the catalytic rate to produce protein P. Since the repressor P r , freely diffuses, the dimerization reaction belongs to Case I. The gene D is spatially fixed and it is repressed by the freely diffusing c 1 , thus this interaction falls under Case II. We assume that the total concentration of D is conserved, so that D T (x) = D(t, x) + c 2 (t, x). The reaction diffusion equations corresponding to (41) are where v Pr (x) (χ r = D r /(L 2 µ)) and v c1 (x) (χ c = D c1 /(L 2 µ)) are the available volume profiles (dimensionless diffusion coefficients) of P r and c 1 , respectively, and from Equation 10 in the Main Text, v c1 (x) = v 2 Pr (x). The boundary conditions corresponding to (42) are Values for dimensionless parameters: The growth rate we used to nondimensionalize the time scales was µ = 0.5 1/hr, consistent with the experiments [10]. The length of the cell is about 3µm and thus L = 1.5µm [8]. The diffusion coefficient of the transcription factor is taken to be D r = D c1 = 0.4µm 2 /s (that of LacI) [13] and thus χ r = χ c = 1280. The transcription factor was assumed to be stale thus γ r = µ. The total concentration of D given byD T was used to nondimensionalize the other concentration variables such that D T = 1. Additional simulation details for Figure 6 in the Main Text: D T (x) ∝ e −20x when DNA near mid cell and D T (x) ∝ e 20(x−1) when DNA at cell poles. The transcription factor production was kept roughly spatially constant α s (x) ∝ e −.001x . The binding and unbinding coefficients were chosen to be a 1 = a 2 = 1000 and d 1 = d 2 = 1000 such that the dissociations constants K d,1 = d1 a1 = K d,2 = d2 a2 = 1.

Approxmating BCF from known parameter values:
From Equations 2 and 10 in the Main Text, we observe that the effective radius of gyration of a dimer complex is √ 2r, where r is the radius of gyration of the individual species. In [14] it was estimated that the radius of gyration for the Tet repressor dimer is 3.1 nm and thus we estimate the radius of gyration of the monomer as r = 3.1/ √ 2. From the expression for θ * given by Equation 26 in the Main Text, Figure 4-B in the Main Text, and r = 3.1/ √ 2, we have that θ * ≈ 0.99 and θ * ≈ 1.01, when the target DNA is near mid-cell Figure K: Infinite dimensional dynamics decay in time independently of binding/unbinding speed Let y(t, x) be as in (42) (represents P r (t, x) or c 1 (t, x)) and let y ⊥ (t, x) = y(t, x) −ȳ(t)v(x) be the a measure of the error in our approximation, where v(x) is the available volume profile of the species. Let ||y(t, x)|| = ( 1 0 (y 2 (t, x))) 1/2 , η 1 = d 1 /µ, and η 2 = d 2 /µ and recall that in the simulations we keep d 1 /a 1 = d 2 /a 2 = 1. We show ||y ⊥ ||/||y||, the relative error for several values of η 1 and η 2 over dimensionless time (with respect to µ) both when the DNA is near mid-cell and the cell poles. The other simulation parameters are identical to those used to generate Figure 6 in the Main Text. For each time point shown, we took the max relative error with respect to the sizes of P r used to generate Figure 6 in the Main Text. The error is high at t = 0 since the initial spatial profiles were chosen to be a constant but note that they quickly decay. Figure L: The relative error between the space averaged PDE model and the reduced ODE model from the data in Figure 6 in the Main Text The relative error in the steady state space averaged protein for the full PDE model (P P ) and the reduced model (P o ) from the data in Figure 6 in the Main Text. The errors are less than 1.7% both when the DNA is localized at mid-cell and at the cell poles. and the cell poles, respectively. Thus, for TetR, the binding strength between the repressor and the DNA varies by about 1% with respect to a well-mixed model in this parameter range. In [15] it was estimated that the radius of gyration for the Lac repressor tetramer is r = 5.3 nm. Assuming that the tetramer is made up of two dimers, then the radius of gyration of each individual dimer is given by r = 5.3/ √ 2. From the expression for θ * given by Equation 26 in the Main Text, Figure 4-B in the Main Text, and r = 5.3/ √ 2, we have that θ * ≈ 0.97 and θ * ≈ 1.03, when the target DNA is near mid-cell and the cell poles, respectively. Thus, for the Lac repressor the binding strength between the transcription factor and the DNA varies by about 3% with respect to a well-mixed model in this parameter range.
While we could not find an exact value for the radius of gyration of the dcas9-gRNA complex, in [16] it was shown that the size of the Cas9-gRNA complex is roughly 10 nm. If we assume this value to be the radius of gyration, from Figure 4-B in the Main Text, we have that the BCF for this complex is θ * = 0.9 and θ * = 1.1 when the target DNA is near mid-cell and the cell poles, respectively. Thus, the binding strength between the Cas9-gRNA complex and the DNA varies s by about 10% with respect to a well-mixed model in this parameter range.

Oscillator
Now we consider the repressor activator clock genetic circuit [17]. This circuit produces sustained oscillations if tuned within an appropriate parameter range [18,19]. The circuit consists of two proteins P a and P r . Protein P a , is an activator which dimerizes to form P a,2 and then binds to its own gene D a to form complex c a,1 to initiate transcription. The dimer P a,2 also finds to the gene D r , which transcribes P r to form complex c a,2 and initiates transcription. Protein P r , dimerizes to form P r,2 and then represses P a by binding to D a to form complex c r . The biochemical equations corresponding to this circuit are: where a i (d i ) for i = 1, . . . , 5 are association (dissociation) rate constants, γ a (γ r ) is the degradation rate of P a (P r ) κ 1 (κ 2 ) is the basal rate at which gene D a (D r ) is transcribed, and κ 3 (κ 4 ) is the rate at which the DNA-transcription-factor complexes are transcribed for D a (D r ). We assume that the total concentration of D a is conserved, so that D a,T (x) = D a (t, x) + c a,1 (t, x) + c r (t, x). Similarly, we assume that the total concentration of D r is conserved, so that D r,T (x) = D r (t, x) + c a,2 (t, x). The spatiotemporal dynamics describing (43) are given by where v Pa (x), v Pr (x), v Pa,2 (x), and v Pr,2 (x) are the available volume profiles of P a , P r , P a,2 , and P r,2 , respectively, χ a = D a /(L 2 µ) is the dimensionless diffusion coefficient of P a and P a,2 , χ r = D r /(L 2 µ) is the dimensionless diffusion coefficient of P r and P r,2 . From Equation 10 in the Main Text, v Pa,2 (x) = v 2 Pa (x) and v Pr,2 (x) = v 2 Pr (x). The boundary conditions corresponding to (44) are Parameters for Figure 7 in the Main Text: The growth rate we used to nondimensionalize the time scales was µ = 0.5 1/hr, consistent with the experiments [10]. The length of the cell is about 3µm and thus L = 1.5µm [8]. The diffusion coefficient of the transcription factor is taken to be D a = D r = 0.4µm 2 /s (that of LacI) [13] and thus χ a = χ r = 1280. The following dimensionless parameters were chosen such that the well-mixed model displayed sustained oscillations: a 1 = 220, d 1 = 1000, a 2 = 1000, d 2 = 1000, a 3 = 1000, d 3 = 1000, a 4 = 1000, d 4 = 1000, a 5 = 1000, d 5 = 1000, κ 3 = 250, κ 1 = .04, κ 4 = 30, κ 2 = .004, γ a = 1, γ r = 0.5. Furthermore, we choose d i and a i for i = 1, . . . 5 large, to demonstrate our results hold even for large binding and unbinding rates. The total concentration of D a which is the same as D r since we assume they are on the same plasmid, is given byD T and it was used to nondimensionalize the other concentration variables such thatD T = 1. The total DNA spatial profile was chosen as D T (x) ∝ e 50(x−1) to model DNA at cell poles.

Numerical Method
In general, a closed form solution to nonlinear PDEs appearing in this work (e.g., Equation 4 in the Main Text) is not available. Therefore, we rely on numerical solutions to directly integrate the PDE. In particular, we utilize a finite difference method that is widely used to simulate PDEs [20]. For a general diffusing species with its concentration given by y(t, x), and available volume profile v(x), we make the following coordinate transformation u(t, x) := y(t, x)/v(x). In this coordinate, the boundary conditions (if applicable) are Neumann [21]: ∂u ∂x x=0,1 = 0 and thus are simpler to implement. Furthermore, with this transformation we observe less stiffness in the numerical simulations. We discretized the spatial domain into N +1 equidistant points such that ∆ = 1/N . Using a second order finite difference method, we approximate the derivatives as: Figure M: The activator is excluded from the chromosome as its size increases The steady state spatial profiles normalized by the average values for P a and P r , that is P a (∞, x)/P a (∞) and P r (∞, x)/P a (∞) for the results of Figure 7 in the Main Text. As the size of P a increases it is excluded from the chromosome. The repressor remains homogeneously distributed throughout the cell. The parameter values and simulation details are identical to those of Figure 7 in the Main Text. Figure N: The relative error between the space averaged PDE model and the reduced ODE model from the data in Figure 7 in the Main Text The relative error in the steady state space averaged activator protein concentration for the full PDE model (P a,P ) and the reduced model (P a,o ) from the data in Figure 7 in the Main Text. Note that for the case when P a is small, large relative errors occur near whenP a (t) reaches a minimum during each period of oscillation. Otherwise all errors are less than 1%.
that appear in the flux terms (in the case where species freely diffuses). The new boundary conditions give rise to the following constraints: This discretization leads to a system of N + 1 ODEs. The resulting set of ODEs are then simulated with MATLAB, using the numerical ODE solvers ode23s. To calculate the space averaged concentrations, we implement a second oder trapezoidal integration scheme [21] given by: For all simulations in this paper N = 200. Now we demonstrate the convergence rate of our numerical method. For the simulation in Figure 5-A in the Main Text when the DNA is localize at the cell-poles and r s /r * = 1, we varied the number of spatial nodes used to discretized the spatial domain to demonstrate the convergence rate of our numerical scheme. Let N be the number of points used to discretize the spatial domain and letm N (t) as the space averaged mRNA concentration for that given discretization. We considered when N = 1024 to be the true solution m 1024 (t) and thus define the following relative error where we took its max value over the time interval of the simulation t ∈ [0, T ] where T = 100. The results from this numerical experiment are shown in Figure O. The convergence rate of our numerical scheme is O(N −2 ) as expected for a second order finite difference numerical scheme.
Figure O: Convergence rate of numerical scheme used to simulate PDEs in this work For the simulation in Figure 5-A in the Main Text when the DNA is localize at the cell-poles and r s /r * = 1, we varied the number of spatial nodes used to discretize the spatial domain to demonstrate the convergence rate of our numerical scheme. Let N be the number of points used to discretize the spatial domain and letm N (t) as the space average mRNA concentration for that given discretization. We considered when N = 1024 to be the true solution m 1024 (t) and the relative error e ∞ is given by (45). The relative error is given by the red markers and the blue dashed lines serve as references for O(N −2 ) convergence rates.
7 Estimating r * from Concentration Profiles and Estimating the BFC Estimate r * : As discussed in Remark 1 in the Main Text, we expect the concentration profile of a freely diffusing species to mirror that of the normalized available volume profile. That is, for a freely diffusing species y, with concentration y(t, x), and available volume profile v(x), we expect Suppose the radius of gyration of y denoted by r, is known and as discussed in the Main Text, we have that v(x) = e −(r/r * ) 2ρ (x) . Approximatingρ(x) as a step function as in (31), we have that v(1) ≈ 1 and v(0) = e − (r/r * ) 2 2(1−∆x) , where ∆x is the distance between the end of the chromosome and the cell poles (see Figure 1 in the Main Text). Let y in and y out denote the the average concentration inside and outside the nucleoid, respectively, which are given by Let ψ y = y out /y in and from (46), we have that then r * can be estimated assuming one knows ∆x, that is, how far the dense nucleoid region extends beyond mid-cell and the average concentration of a species inside and outside nucleoid region. A similar calculation was done in [8], using the fact that the free ribosome concentration is 10% higher at the cell poles than mid-cell. Estimate the BCF: The BCF provides a measure to determine the extent to which spatial effects modulate the biomolecular dynamics. Therefore, an experimental method to estimate the BCF is desirable. We propose a method that only requires knowing ∆x and the concentration of freely diffusing species inside and outside the nucleoid . Suppose that for Case I and Case III, the concentration of E i is measured inside and outside the nucleoid and denoted by E in i and E out i , respectively. Similarly, for Case I and Case II, we assume that S in and S out is measured. If a fluorescence imaging method is used to measure these quantities (as in [9]), then we emphasize that the free E i and S must be measured, not when they are in complex form (c i ). Let ψ i = E out i /E in i and ψ s = S out /S in . By (47), ψ i = v 0,Ei and ψ s = v 0,S , where v 0,Ei and v 0,S are as in (32). Thus, using (33) we can estimate the BCF for Cases I-III by • Case I: • Case II and Case III: As r e,i /r * → ∞ (r s /r * → ∞) we have that ψ i → 0 (ψ s → 0), thus ψ i and ψ s are a measure of the excluded volume effects on E i and S, respectively. Physically, this is expected because when E i is severely expelled from the nucleoid by available volume effects, we have that E out i E in i and similarly for S.

Experimental Setups to Verify the Role of Spatial Effects Predicted by Model
A potential experiment to test our hypothesis that genes near the poles our transcribed more effectively than gene near mid-cell, is to measure the rate of transcription (via Quantitative PCR) of a gene under the control of the T7 promoter. This promoter is solely transcribed by the T7 RNAP which specifically targets the promoter, thus this system can be considered orthogonal to the endogenous transcription machinery [22]. By appending random base pairs (BPs) to the sequence of T7 RNAP that do not affect its functionally, we can control its size and thus the extent of excluded volume effects. We can then measure the transcription rate of the gene regulate by the T7 promoter when it is localized in the cell-poles and mid-cell. The results of this experiment should look similar to Figure 5-A in the Main Text, as the size of the T7 RNAP increases the mid-cell (pole) gene has lower (higher) transcription rate. For the mid-cell localized gene, this can be repeated in parts of the chromosome which are known to be dense to amplify these effects.
To experientially validate our analytical prediction that protein steady state levels will increase with mRNA size, we propose expressing a florescence protein from a plasmid with an appended sequence of base pairs added downstream of the stop codon. The appended sequence should have a low affinity to recruit ribosome such that the amount of ribosomes sequestered by the mRNA are the same as without the appended sequence. Assuming this appended sequence does not affect the lifetime of the mRNA, then it should yield the same functional protein which can be used to quantify the mRNA excluded volume effects. This appended sequence of base pairs will allow us to control the size of the mRNA without increasing its ribosome usage. From our theory, for longer appended sequences, more protein expression is expected.
To validate the hypothesis that a transcriptional repressor regulates genes near the poles more effectively than gene near mid-cell, we propose a genetic circuit on a plasmid expressing a repressor that targets a gene expressing a florescence protein. The transcription factor chosen should be large enough or dimerize to have considerate excluded volume effects. The target DNA expressing protein should be placed on several axial locations in the cell (under the same promoter) achieved by using backbones with different localization profiles and/or different chromosomal integration sites. We should observe that the effective disassociation constant of the repression curve increases as the target genes location is closer to the mid-cell. The disassociation constant is proportional to the amount of repressor necessary to cause the genes expression to decrease by half.

9
Cell Division: Time Varying Cell Length and Chromosome Profile As the cell divides it partitions molecular species amongst daughter cells, this along with changes in the cell length cause dilution effects on intracellular concentrations. Furthermore, early in the cell division cycle, the chromosome density is highest mid-cell, but as the cell divides the peak chromosome density tends towards the cell-poles [9,8] (to distribute genes evenly among daughter cells). From the results in the Main Text, we expect this temporal changes in the chromosome density will effect the BCF since species are repelled away from regions with high chromosome density via excluded volume effects. In this section, we provide the modeling framework to account for dilution effects and temporal fluctuations in the chromosome density.
Dilution effects on spaced average concentrations: Here we demonstrate how cell division and a time varying cell length effects space average concentrations. LetN p (t) be the total molecular count in a cell population of a molecule of interest (i.e., ribosomes) as the cell expands and divides. To model dilution from cell division, assume thatN p (t) is identically distribute among N cells (t) number of cells such that N cells (t) = N cells (0)e µt , where µ is the cell growth rate and each cell has a volume given by V c (t) = 2πR 2 L(t), where R is the cell radius and 2L(t) is the cell length. The total population volume is then given by V (t) = N c (t)V c (t) and lettingc(t) be the number of molecules per total volume (concentration), this quantity is given byc(t) =N p (t) V (t) , note that this is identical to the concentration per cell volume (since we assume all cells in the population have identical averaged concentrations). This implies thaṫ Dilution effects on local concentrations: Let N (t, x) be the number of molecules per unit length of a cell such that that . The temporal evolution of N (t, x) in the presence of dilution and a moving boundary, which introduces an advective term [23] (to account for the extra diffusion as the cell length varies), is given by with boundary conditions where u(t, x) = x L(t)L (t) is the velocity of a material point induced by the increase in cell length. Notice that our current finite difference method with a stationary mesh cannot be applied directly to (49), thus we propose the following spatial coordinate transformation y(t, x) = x L(t) (and to be consistent with the nondimensionalization from the Main Text), which renders a stationary domain. Let c(t, y) := N (t, yL(t)) and thus dN dt = dc dt = ∂c ∂t + ∂c ∂y ∂y ∂t = ∂c ∂t − ∂c ∂y yL L , finally Notice that the effective dilution coefficient is now given by D L 2 (t) , which as expected increases as cell length increases and 2 1 0 c(t, y)dy = 2 L(t) dx =ċ(t), thus the space averaged under this coordinate system provides us the concentration per cell volume. The boundary conditions are Time varying chromosome density: We now model the chromosome density varying in timeρ :=ρ(t, x) as the cell divides. This implies that the available volume profiles will also depend on time since v(x, t) = e − r r * ρ (x,t) , and thus ∂c ∂t whereμ(t) is the effective dilution rate given bỹ The quantities L(t) and v(t, x) will vary with a time scale related to cell growth, for example, let T 1/2 = ln (2) µ be the cell doubling time, then one possibility is for this choice of L(t), the effective dilution rate (50) is graphically shown in Figure P. In [8] it was shown the cell length late in the cell division cycle was 4.4µm (compare to its nominal length 3µm), thus for our simulations we take ∆ L = 0.2.
Figure P: Varying cell length modulates dilution rate The effective dilution rate (50) is given for L(t) given by (51) and ∆ L = 0.2, where t 1/2 is time normalized by the cell doubling time (t 1/2 = t/T 1/2 ).
in time when performing model reduction as in the "Time scale separation" section in the Main Text, thus we expect (similar to the results of the Main Text) wherec(t) = 1 0 c(t, y)dy andv(t, y) = v(t, y)/ 1 0 v(t, y)dy. So all of our previous analysis still holds except that the BCF will vary slowly (slow with respect to the time scale of diffusion) as the cell divides. Example: We verify via simulations the prediction that (53) holds and that the BCF can be treated as a slowly (with respect to diffusion) varying parameter. Consider the simple bimolecular reaction: with dynamics given by ] + α e (y) − aE(t, y)S(t, y) + dc(t, y) −μ(t)E(t, y), where α e (y) and α s (y) are the production rates of E and S, respectively. The space averaged dynamics (Ē(t) = to show the effects of having time varying dilutionμ(t) and available volume profiles on the expression level. Figure R shows how the space averaged concentrationĒ(t) is modulated by the time varying dilutionμ(t).
We observe that the concentration reaches a periodic steady state centered at unity where the oscillations have a period that coincides with the doubling time. Furthermore, in Figure S, we verify that wherev e (t, y) = v e (t, y)/ 1 0 v e (t, y)dy as expected (since diffusion much faster than dilution D e /(L 0 µ) 1). Thus, even with a time varying v e (t, y), the enzyme will be expelled from the chromosome to areas of higher available volume. Next, we demonstrate how the binding dynamics are effected by having a time varying available volume profile, therefore a = 0 in (55). Similar to "Case 1" in the Main Text, we consider the case when D e , D s , D c = 0 (all species freely diffuse), where we expect the BCF to be approximated by θ * (t) =   This is verified in Figure T, where θ(t) given by (57) and θ * (t) given by (59) are shown after three doubling times and are shown to be in good agreement. The BCF varies periodically in time (with the period consistent with the doubling time) and oscillates near a nominal value of 1.5 with amplitude 0.04. Next we look at the the case when S and c are spatially fixed (D s = D c = 0) and localized near y * , which is similar to "Case 2" in the Main Text. For this scenario we expect the BCF to be approximated by θ * (t) =v s (t, y * ).
The results are shown in Figure U when y * = 0 and y * = 1 after three doubling times, there is good agreement between the BCF and its approximation. When y * = 0 (S localized near mid-cell), the BCF is less than unity and oscillates near a nominal value of 0.55 with amplitude 0.3. When y * = 1 (S localized near the cell poles), the BCF is greater than unity and oscillates near a nominal value of 2 with amplitude 0.1. These results suggest that the BCF for a species localized near mid-cell will vary significantly as the cell density varies during cell division. This is expected because as shown in Figure 1 in the Main Text, the chromosome density is initially high near mid-cell but decreases by half as the cell divides, thus no longer excluding from that region.

Exclusion Effects from Plasmid DNA Density
The genome of E. coli MG1655 has 4.6 Mbp [24]. Comparatively, a single plasmid can have .01 Mbp a copy number as high as 500-700 (e.g. pUC19). Therefore, the total plasmid and chromosome basepair count may be comparable in applications with high copy number plasmids. In these applications, it may be necessary to account how plasmid DNA repels freely diffusing species and "excludes" them. To do so we modify our model of the DNA densityρ(x) as shown in Figure V to account for plasmid DNA. For the DNA density profiles from Figure V, we calculate the approximate BCF θ * (Equation 5 in the Main Text), these are shown in Figure W. For Case 1 where the reactants freely diffuse. We observe that the BCF decreases as the plasmid DNA density increases (as shown in Figure V). This occurs because as the plasmid DNA increases the overall density profile becomes more uniform. Note that when the plasmid DNA is sufficiently high to render an almost uniform DNA density profile, the BCF is unity as expected. For Case 2, where one reactant freely diffuses and the other is fixed at x * . As the plasmid density increases (as shown in Figure V) we observe that the BCF decreases at the cell poles (as expected since species are excluded from dense plasmid DNA mesh) and increases at region near x * ≈ 0.65 (where there is minimal overlap between chromosome and plasmid DNA). When the plasmid DNA density is similar to that of the chromosome rending a uniform DNA distribution, we observe that the BCF is unity everywhere (as expected since there are no exclusion effects).
(a) y * = 0 (b) y * = 1 Figure U: The BCF for the case when S is spatially fixed at y * The binding correction factor θ(t) (57) and its approximation θ * (60) over two cell division cycles. The BCF oscillates with a period consistent with the doubling time. When y * = 0 (S localized near mid-cell), the BCF is less than unity and oscillates near a nominal value of 0.55 with amplitude 0.3. When y * = 1 (S localized near the cell poles), the BCF is greater than unity and oscillates near a nominal value of 2 with amplitude 0.1. This is shown for a molecule localized near mid-cell y * = 0 and near the cell poles y * = 1. The simulation parameters are r e /r * = r s /r * = 2, r c /r * = 2 √ 2, µ = 1, D e /L 0 = D s /L 0 = D c /L 0 = 13 × 10 3 α e (y) = α s (y) = 1, ∆ l = 0.2, d = 100, a = 100. where x s ∈ [1/2, 3/2] is the parameter we varied to get difference plasmid densities (x s = 3/2 lowest plasmid density and x s = 1/2 highest plasmid density).

Codes
A copy of the codes used in this study can be found at: https://drive.google.com/drive/folders/11tJ9rn8QLs-G2xOsOhpXnOSqPJvI804y?usp=sharing The approximate BCF θ * when plasmid DNA is accounted for (a) Case 1 where the reactants freely diffuse. We observe that the BCF decreases as the plasmid DNA density increases (as shown in Figure V). This occurs because as the plasmid DNA increases the overall density profile becomes more uniform. Note that when the plasmid DNA is sufficiently high to render an almost uniform DNA density profile, the BCF is unity as expected. (b) Case 2 where one reactant freely diffuses and the other is fixed at x * . As the plasmid density increases (as shown in Figure V) we observe that the BCF decreases at the cell poles (as expected since species are excluded from dense plasmid DNA mesh) and increases at region near x * ≈ 0.65 (where there is minimal overlap between chromosome and plasmid DNA). When the plasmid DNA density is similar to that of the chromosome rending a uniform DNA distribution, we observe that the BCF is unity everywhere (as expected since there are no exclusion effects). The simulation parameters are r/r * = 1 andρ(x) as in Figure V.