Optimal control of stochastic reaction networks with entropic control cost
and emergence of mode-switching strategies

Let us consider a population of particles that evolve through stochastic events, i.e., the occurrence of reactions. Each particle is assigned one of discrete type set $X$. The number of particles of type $x\in X$ at time $t$ is denoted as $n_{x}(t)\in\mathbb{Z}_{\geq 0}$. The number distribution $n(t)\in\mathbb{Z}_{\geq 0}^{|X|}=\mathcal{N}$ completely characterize the state of the system at time $t$. The number distribution changes when a reaction occurs. There are $|R|$ kinds of reactions, where $R$ is the set of reactions. Once a reaction $r\in R$ occurs, $s^{-}_{r,x}\in\mathbb{Z}_{\geq 0}$ particles are consumed and $s^{+}_{r,x}\in\mathbb{Z}_{\geq 0}$ particles are produced. The net change in the number of particles $s_{r,x}:=s_{r,x}^{+}-s_{r,x}^{-}\in\mathbb{Z}$ is called stoichiometric coefficient and $s_{r}=(s_{r,x})_{x\in X}\in\mathbb{Z}^{|X|}$ is stoichiometric vector for reaction $r$. Thus when the reaction $r$ occurs at time $t$, the number distribution changes from $n(t^{-})$ to

The timing of reaction events is random, which follow an inhomogeneous Poisson process with rate $\lambda_{r}(t)\geq 0$ for reaction $r$ at time $t$. It could vary over time and depends on the current number distribution $n(t)$. We write the rate of reaction $r\in R$ as the product of the reaction rate coefficient $k_{r}(t)\geq 0$ and the propensity function $h_{r}(n)\geq 0$:

When we assume the law of mass action kinetics, $k_{r}(t)$ is a kinetic constant, and the function $h_{r}(n)$ is given by $h_{r}(n)=N^{1-\sum_{x\in X}s^{-}_{r,x}}\prod_{x\in X}\frac{n_{x}!}{(n_{x}-s^{-}{r,x})!}$, where $N$ is a constant parameter that describe the system size, e.g., the volume or total number of the particles in the system. We do not assume the mass action kinetics in the following so that our results become applicable to any propensity function. Instead, we presume only the situations that make the process well-defined up to time $T$ or up to exit time $T_{\mathrm{exit}}$. For example, $h_{r}(n)$ should be $0$ when the reaction pushes the state out of the domain, i.e., $n+s_{r}\not\in\mathcal{N}$.

The stochastic reaction system defined in this way can be characterized by either the counting process representation or the Markov process representation. The counting process representation [20] is useful for characterizing the stochastic process $(n(t))_{0\leq t\leq T}$ and simplifies the derivation of KL control cost. Let $\xi_{r}(t)\in\mathbb{Z}_{\geq 0}$ denote the number of occurrences of reaction $r$ from time $0$ to $t$. The number distribution $n(t)$ is a linear transformation of the reaction count $\xi(t)$, i.e., for any $0\leq t\leq T$,

The probability law of $\xi(\cdot)$ is given by the Poisson processes. Using the time-change of independent unit Poisson processes $Y_{r}(t)$, the reaction count can be written as $\xi_{r}(t)=Y_{r}\left(\int_{0}^{t}k_{r}(\tau)h_{r}(n(\tau))d\tau\right)$.

The stochastic process $n(\cdot)$ is a Markov process on the nonnegative integer lattice $\mathcal{N}$. The transition rate from $n\in\mathcal{N}$ to $n^{\prime}\in\mathcal{N}$ is given by

where $\delta_{n,n^{\prime}}$ is the Kronecker delta. Let $P_{t}(n^{\prime}):=\mathbb{E}\left[\delta_{n(t),n^{\prime}}|n(0)=n_{0}\right]$ denote the probability of being the state $n^{\prime}\in\mathcal{N}$ at time $t$ given the initial state $n_{0}\in\mathcal{N}$. Then, $P_{t}(n^{\prime})$ satisfies the Kolmogorov forward equation (chemical master equation) for any $0<t<T$ and $n\in\mathcal{N}$:

We assume that the controller can observe the current state $n(t)$ and adjust the reaction rate coefficients $k_{r}(t)$ to any desired value at any time while the function $h_{r}(n)$ remains unchanged, i.e., the controller can modulate only the speeds of reactions. We would like to find the optimal reaction rate coefficients $(k(t))_{0\leq t<T}$ which drive the population to take a desirable trajectory $(n(t))_{0\leq t<T}$ while minimizing the cost of modulating the reaction rate coefficients. Let us define the utility function $\mathcal{U}(n(\cdot)):=U_{T}(n(T))+\int_{0}^{T}U_{\tau}(n(\tau))d\tau$ as the sum of the terminal utility function $U_{T}(n)\in\mathbb{R}$ and the time integral of instantaneous utility function $U_{\tau}(n)\in\mathbb{R}$. The control cost function is also defined as $\mathcal{C}(n(\cdot),k(\cdot)):=\int_{0}^{T}C(n(\tau),k(\tau))d\tau$. Then, we consider the following optimization problem:

or equivalently, the following unconstrained optimization:

where the expectation $\mathbb{E}_{k(\cdot)}$ is taken over trajectories $n(\cdot)$ generated under the designated reaction rate coefficients $k(\cdot)$ and initial condition $n(0)=n_{0}$. A scalar $\beta$ is the Lagrange multiplier or a parameter to adjust the importance of the utility relative to the control cost.

We will calculate the optimal transition rate $(k^{\dagger}(t))_{0\leq t<T}$ that attains the optimum of $\mathcal{E}(n_{0},u)$ and $V(n_{0},\beta)$. The standard protocol for solving the problem is to consider the value function $V_{t}(n,\beta)$, which is defined as the maximum of the expectation in Eq. (7) from $t$ to $T$ under the condition $n(t)=n$:

where $\mathcal{U}_{t}(n(\cdot))$ and $\mathcal{C}_{t}(n(\cdot),k(\cdot))$ are the utility and control cost functions from time $t$ to $T$, e.g., $\mathcal{U}_{t}(n(\cdot)):=U_{T}(n(T))+\int_{t}^{T}U_{\tau}(n(\tau))d\tau$. The value function is equal to the terminal utility function at time $T$, $V_{T}(n,\beta)=\beta U_{T}(n)$, and also provides the optimum for the original problem at time $0$, $V_{0}(n,\beta)=V(n,\beta)$. Thus, the transition rate attaining Eq. (8) is identical to the optimal transition rates $k^{\dagger}$ in the original control problem.

If $V(n,\beta)$ were analytically or numerically obtained for all $n$ and $\beta$, $\mathcal{E}(n,u)$ in Eq. (6) is derived as the Legendre transformation of $V(n,\beta)$: (see Sec. A in the Supplementary Material ¹¹1See Supplemental Material at for the derivation of equations and detailed numerical results, which includes Refs. [59, 60, 61, 62, 63, 64]. for the derivation)

Furthermore, the expected utility under the optimal control can also be derived from the value function (see Sec. A of the Supplementary Material [27] for the derivation)

To obtain the value function, we usually leverage the Hamilton–Jacobi–Bellman (HJB) equation, a differential equation that $V_{t}(n,\beta)$ satisfies. However, solving the HJB equation is generally intractable both analytically and even numerically because it is a nonlinear differential equation on a possibly infinite domain. This difficulty is the major obstacle to optimal control in applications. To address the difficulty, we are conventionally forced to restrict or approximate the original problem to a linear quadratic Gaussian (LQG) problem, in which only the linear dynamics on continuous Euclidean space with Gaussian noise with a quadratic control cost function can be considered. However, all these restrictions of the LQG problem conflict with the essential properties of stochastic RN control, i.e., the nonlinearity of $h_{r}(n)$, the discreteness and nonnegativity of state $n\in\mathbb{Z}_{\geq 0}^{|X|}$, the nonnegativity of the control parameter $k_{r}\geq 0$, and Poissonian nature of stochasticity. Several studies attempted to overcome a part of these difficulties [28, 29, 30]. Nonetheless, optimal control was not practical for biological problems described by RN.

In this work, we clarify that the difficulty can be resolved by adopting the following cost function

where $c_{KL}(a,a^{0})$ is the generalized Kullback–Leibler (KL) divergence for positive scalars $a,a^{0}>0$:

The value at $a=0$ is defined by the limit $c_{KL}(0,a^{0}):=\lim_{a\rightarrow 0}c_{KL}(a,a^{0})=a^{0}$.

The KL control cost function in Eq. (11) measures the deviation of the controlled reaction rate coefficient $k$ from its uncontrolled reaction rate coefficient $k^{0}\in\mathbb{R}^{|R|}_{>0}$. Because the control cost function $C(n,k)$ is minimized at $k=k^{0}$, the optimal reaction rate coefficient $k^{\dagger}$ is equal to $k^{0}$ when $\beta=0$, i.e., no preference on the state trajectory. If $\beta\neq 0$, the optimal reaction coefficient for $r\in R$ at time $t$ is given by

where $\overline{\nabla}_{s_{r}}f(n):=f(n+s_{r})-f(n)$ is the discrete gradient (see Sec. A of the Supplementary Material [27] for the derivation). As the slope of $c_{KL}$ goes to infinity, $\frac{\partial}{\partial k_{r}}c_{KL}(k_{r},k_{r}^{0})\rightarrow-\infty$ for $k_{r}\rightarrow 0$, the nonnegativity constraint over $k$ is automatically satisfied. In fact, the exponential function in Eq. (13) ensures that the optimal coefficient $k^{\dagger}_{r}$ is non-negative.

Note that the KL control cost in Eq. (11) is weighted by the propensity function $h_{r}(n)$ for each reaction $r\in R$, which makes the control cost function depend on the current state $n$. It means that even if the reaction rate coefficient is the same, the control cost for a reaction $r$ is higher when the reaction occurs more frequently. A similar form of the control cost function has been used in quantifying the biological cost of cellular reactions [31].

For the KL cost function, the HJB equation becomes

which is yet a nonlinear differential equation. However, by the Cole–Hopf transformation (logarithmic transformation), $Z_{t}(n,\beta):=\exp(V_{t}(n,\beta))$, the HJB equation is linearized as

with the terminal condition $Z_{T}(n,\beta)=\exp(\beta U_{T}(n))$.

Although it is still an ordinary differential equation on a possibly infinite domain, linearity allows us to calculate the value function efficiently. Note that the second term on the right-hand side is the backward operator (generator) $\mathcal{L}_{k}$, which is the adjoint of the forward operator $\mathcal{L}^{\dagger}_{k}$ in the chemical master Eq. (II.1). Moreover, the special form of Eq. (15) allows us to have the following probabilistic representation for $Z_{t}(n,\beta)$:

which is called the Feynmann-Kac formula, cf., Appendix 1. Prop. 7.1 of [32]. Thus, the value function is given by

This representation enables us to compute the value function by evaluating the expectation of the utility function with respect to the uncontrolled reaction rate coefficient $k^{0}$. For simple reaction systems, we could obtain the analytical expression of the value function as demonstrated in the following sections. Even if it is not possible, efficient Monte Carlo sampling techniques can be used to estimate the expectation. It is worth noting that this representation of the value function can be evaluated in a time-forward manner, whereas the usual optimal control problem requires the time-backward calculation of the HJB equation due to the dynamic programming principle.

The linearization of the HJB equation and efficient computation of the optimal control problems for certain control cost functions was reported by Kappen [24] for diffusion processes and by Todorov [26] for discrete-time Markov chains. As Theodorou and Todorov [33] discussed, the linearization of the HJB equation is possible if the control cost function is given by KL divergence between path measures. In this case, the optimal control problem is related to an estimation (filtering) problem. Similar mathematical properties had been found in relation to the duality of control and inference [21, 22, 23]. Similar properties for Markov jump processes are identified in recent studies [34, 35, 36] as well as one of the earliest studies by Fleming [21]. In this paper, we develop the optimal control framework for stochastic RN inspired by these previous studies. In Sec. B of the Supplementary Material [27], we elaborate on the path measure perspective and why the KL cost function in Eq. (11) works.

We have hitherto focused on the problem with a finite and fixed terminal time $T$. One can extend our approach to the problem with infinite time length $T\rightarrow\infty$. Since a stochastic reaction system often has absorbing states $\mathcal{N}_{\mathrm{abs}}\subset\mathcal{N}$, the stochastic process inevitably reaches one of these absorbing states after some time $T_{\mathrm{exit}}:=\min\{t\geq 0|n(t)\in\mathcal{N}_{\mathrm{abs}}\}$. For example, in the birth-death processes

the extinction state $n_{A}=0$ is absorbing. The extinction events are particularly important in biological problems because the extinction of some spieces in a chemical, organismal, or human population could alter the dynamics of the system qualitatively. Thus, controlling species into either survival or extinction has many applications, as mentioned in the Introduction. Moreover, biological control problems may not always have a prescribed end time $T$ because the goal of control is usually to achieve something rather than to do something until $T$. Thus, the first exist control is more essential than fixed-time control. Nevertheless, it has been less practiced because the first exist control is more involved technically.

We can consider the optimal control problem with such random terminal time $T_{\mathrm{exit}}$ and KL cost. The utility function is replaced with the sum of the terminal utility function $U_{\mathrm{exit}}(n)$ defined on the absorbing states $\mathcal{N}_{\mathrm{abs}}$, and the instantaneous utility function $U(n)$ defined on the non-absorbing states $\mathcal{N}_{\mathrm{non}}:=\mathcal{N}\backslash\mathcal{N}_{\mathrm{abs}}$, i.e.,

Then, the value function becomes time-independent $V(n,\beta)=V_{t}(n,\beta)$. If we assume the KL control cost, and the Cole–Hopf transformation $Z(n,\beta):=\exp(\beta V(n,\beta))$ yields

for $n\in\mathcal{N}_{\mathrm{non}}$, and $Z(n,\beta)=\exp\left(\beta U_{\mathrm{exit}}(n,\beta)\right)$ for $n\in\mathcal{N}_{\mathrm{abs}}$. Thus, the first exit optimal control problem is reduced to a linear algebraic problem. Similar results have been obtained for discrete-time Markov chains [26]. The probabilistic representations in Eqs. (16), (17) are also applicable.

If there are no absorbing states, the terminal time and the time accumulated objective function could diverge. For such cases, the time-averaged formulation is useful:

Via the Cole–Hopf transformation, the optimal solution could be cast into an eigenvalue problem (see Sec. C of the Supplementary Material [27]).

We demonstrate the effectiveness of our method by applying it to several control problems. The first example is the jump processes on discrete states $X$. We consider $N$ particles walking on a directed graph $(X,R)$, where the set of directed edges $R$ represents the allowed transitions between types $X$. For any $x,y\in X$, $(x,y)\in R$ implies that the particle can jump from $x$ to $y$. The stoichiometric coefficients are given by $s_{(x,y),z}=\delta_{y,z}-\delta_{x,z}$ and the kinetics is the mass action type $h_{(x,y)}(n)=n_{x}$ for $x,y,z\in X$.

For $X=\mathbb{Z}$ and $N=1$, the process is reduced to a simple one-dimensional random walk by a single walker. Such a process has been used as a simple model of intracellular transport of macromolecules along intracellular filaments [37]. Molecular motors such as dynein and kinesin consume energy and move in one direction. Since the step length is fixed, the position of a molecular motor on a filament can be modeled as a 1-dimensional discrete grid. For this case, we can have analytical solutions for some optimal control problems owing to the sound properties of KL control.

Consider the situation where a particle is required to reach a goal $x^{*}\in\mathbb{Z}$ as soon as possible. The position of the particle at time $t$ is denoted as $\hat{x}(n(t))\in X$. We can formulate the situation as the following minimum exit time problem:

where the absorbing states are given by $\mathcal{N}_{\mathrm{abs}}=\{n\in\mathcal{N}|\hat{x}(n)=x^{*}\}$. Setting $U_{\mathrm{exit}}(n)=0$ and $U(n)=-1$, we have an equivalent maximization problem

Equation (17) shows that the value function $V(n,\beta)=V(\hat{x}(n),\beta)$ is equal to the cumulant-generating function of the exit time with parameter $-\beta$ and given by

Assuming the symmetric uncontrolled transition rates $k^{0}_{r}=\kappa$ for all $r\in R$ and using the analytical expression of the cumulant-generating function [38], we obtain the value function analytically:

where a scalar $\gamma(\beta)\geq 0$ is defined as

The conjugate value function can be calculated for $x^{*}\neq x$ and $u<0$ as follows

where $\Delta:=x^{*}-x$. Since the value function is linear in $|x^{*}-x|$, we have the piecewise constant optimal transition rates:

When $x<x^{*}$, the transition rate $k^{\dagger}_{(x,x+1)}$ in the positive direction is higher than the uncontrolled one $\kappa$, while the transition rate $k^{\dagger}_{(x,x-1)}$ in the negative direction is lower than $\kappa$. We can also obtain the analytic solution for the standard control problem to maximize the average speed of the walker (see Sec. D of the Supplementary Material [27]).

Next, we consider the case where $N=2$ and $X=\mathbb{Z}$. Let us consider the minimum time problem with exclusion interaction, i.e., all the particles should reach the goal $x^{*}\in X$ as soon as possible while they are not allowed to occupy the same site. This is a model of two molecular motors moving on the same filament on which they cannot overtake (Fig. 1 (a)). Non-colliding random walks and diffusion processes have been studied persistently [39, 40, 41].

Let us denote the position of the left particle as $\hat{x}_{1}(n)$ and the right particle as $\hat{x}_{2}(n)$, i.e., $\hat{x}_{1}(n)<\hat{x}_{2}(n)$. Due to the exclusive interaction, the first particle cannot overtake the second. Thus, the particles are identifiable for all $t\geq 0$. The absorbing states are given by

Among the abosrbing states above, we are interested in the single state $n^{*}\in\mathcal{N}_{\mathrm{abs}}$ such that the particles reach the goal $\hat{x}_{1}(n^{*})=\hat{x}_{2}(n^{*})=x^{*}$. We can formulate the problem as a first exit problem with the instantaneous utility $U(n)=-1$, and the terminal utility

Assuming that each particle stops when it reaches the goal, the time-integrated instantaneous utility is given by

where $T^{i}_{\mathrm{exit}}:=\inf\{t\geq 0|\hat{x}_{i}(n(t))=x^{*}\}$ is the first exit time of the $i$-th particle ($i=1,2$). Then, the value function $V(n,\beta)=:V(\hat{x}_{1}(n),\hat{x}_{2}(n),\beta)$ can be calculated as

where $V(x_{i},\beta)$ is the value function in Eq. (25) for a single random walker. When the first particle is close to the goal $x_{1}\approx x^{*}$ while the second particle is still away from it, $|x_{2}-x^{*}|\gg 1$, the upper bound become tight and we have

Using the upper bound, the conjugate value function $\mathcal{E}(x_{1},x_{2},u)$ satisfies

where $\mathcal{E}(x_{i},u)$ is the conjugate value function in Eq. (27).

Analytical estimates are compared with numerical solutions in Fig. 2. Numerical solutions are obtained by solving the linear equation in Eq. (20) where the infinite space $\mathbb{Z}$ is truncated to the finite interval $X=\{0,1,\ldots,M-1\}$. The goal position is set to $x^{*}=0$.

When the first particle has already reached the goal $x_{1}=0=x^{*}$, the problem is equivalent to the minimum time problem for the second particle only. Then, the value function satisfies $V(x^{*},x_{2},\beta)=-\gamma(\beta)|x^{*}-x_{2}|$, which is consistent with Eq. (33) (Fig. 2 (b), upper panel). Then, the gradient of $V(x^{*},x_{2},\beta)$ is constant (Fig. 2 (c), upper panel), which results in the constant optimal control.

When the first particle is at $x_{1}=10$, the value $V(10,x_{2},\beta)$ as a function of $x_{2}$ has a gentle slope than $\gamma(\beta)$ to avoid collision with the first particle (Fig. 2 (b), lower panel). Especially when $\beta$ is small, collision avoidance is more important than early arrival. In this case, the value gradient can be positive, and the second particle moves away from the first particle, as seen in Fig. 2 (c) lower panel. On the other hand, when $\beta$ is large, the exit time becomes more important than the collision. Thus, the value gradient is always negative except at the collision point, i.e., at $x_{2}=10$.

In the control of birth and death processes of cells or organisms without immigration, the extinction state $n=0$ is a natural absorbing state. In the context of biodiversity conservation, one has to avoid the extinction. Furthermore, a single species should not dominate the ecosystem. To address this population control problem, let us consider the following birth and death reactions involving two species $A$ and $B$ (Fig. 1 (b)):

Assume that every dying individual is replaced by a duplicated individual of either $A$ or $B$, and that the size $N$ of the population is constant. This is the continuous-time Moran process studied in population genetics [42, 43, 44]. Since birth and death reactions occur simultaneously, we can summarize the four reactions into two:

where the propensity functions are given by $h_{1}(n)=h_{2}(n)=n_{A}n_{B}/N$. We denote the ratio between the rate coefficient by $\gamma:=k_{1}^{0}/k_{2}^{0}$, which means that type $A$ has a selective advantage over type $B$ when $\gamma>1$ and vice versa. The case with $\gamma=1$ is the neutral situation. In any case, the system eventually reaches one of two extinction boundaries $(n_{A},n_{B})=(N,0)$ or $(0,N)$ due to the random drift. We formulate the maintenance of the population coexistence as the maximization of the exit time problem from coexisting states:

The value function, i.e., the cumulant generating function of the extinction time, can be decomposed as

where $K_{0}(n,\beta)$ and $K_{N}(n,\beta)$ are the cumulant-generating function of the first hitting time at $0$ and $N$, respectively. They have an explicit expression [45] using the eigenvalues of submatrices of the transition rate matrix $M:=(-\omega_{n,n^{\prime}})_{n_{A},n_{A}^{\prime}\neq 0,N}\in\mathbb{R}^{(N-1)\times(N-1)}$ as follows:

and

where $b^{0}(a):=k_{1}^{0}h_{1}(a,N-a)$ and $d^{0}(a):=k_{2}^{0}h_{2}(a,N-a)$. In the equations, $x_{i}(n_{A})>0$, $y_{j}(n_{A})>0$, and $\lambda_{k}>0$ are the eigenvalues of bottom-right submatrix of size $(N-n_{A}-1)$, top-left submatrix of size $(n_{A}-1)$, and the full matrix $M$, respectively [45]. Due to the interlacing property [46], the smallest eigenvalue $\lambda_{1}$ of the full matrix is smaller than the smallest eigenvalues, $x_{1}(n)$ and $y_{1}(n)$, of the submatrices. Therefore, the value function diverges $V(n,\beta)\rightarrow+\infty$ as $\beta\rightarrow\lambda_{1}$. The expected extinction time $\mathcal{U}^{\dagger}(n,\beta)$ under the optimal control also diverges as $\beta\rightarrow\lambda_{1}$ because it is the derivative of $V(n,\beta)$ by $\beta$ (Eq. (10)). This critical $\beta_{c}:=\lambda_{1}$ is the same for all initial conditions $n$.

We numerically calculated the value function in two ways: by solving the linear algebraic Eq. (20) and by using the analytical formula in Eqs. (38)–(40). Two solutions exactly match, and the value function $V(n,\beta)$ diverges as $\beta$ approaches the smallest eigenvalue $\beta\rightarrow\beta_{c}=\lambda_{1}$ (Fig. 3 (a)). The cost-exit time tradeoff curves $\mathcal{E}(n,u)$ in Fig. 3 (b) approach lines with slope $\beta_{c}$ in the long exit time regime. Thus, the cost rate per unit time $\mathcal{E}(n,u)/u$ is upper bounded by $\beta_{c}$ as shown in Fig. 3 (c). The result indicates that a finite amount of control cost per time is sufficient to prevent the ultimate extinction on average. As $\gamma=k^{0}_{1}/k^{0}_{2}$ increases, the extinction tends to happen earlier, resulting in the increased control cost rate $\beta_{c}$ for preventing extinction.

As a function of $n_{A}$, the value function $V(n,\beta)$ for each $\beta$ and $\gamma$ has a single peak, which we designate by $n_{A}^{*}$ (Fig. 3 (d)). The optimal control steers the system to stay around $n_{A}^{*}$. As $\beta$ increases, the value gradient becomes steep around $n_{A}^{*}$ (Fig. 3 (e)) so that the controlled trajectories do not leave there even after a long time (Fig. 3 (f)).

When $\beta$ is close to $\beta_{c}$, the width of the peak around $n_{A}^{*}$ decreases as $\gamma$ increases, and there emerges a region with a shallow gradient $\overline{\nabla}_{s_{1}}V(n,\beta)\approx 0$. In this region, the optimal control $k^{\dagger}_{r}(n)=k^{0}_{r}\exp(\overline{\nabla}_{s_{r}}V(n,\beta))$ become close to the uncontrolled rate $k^{0}_{r}$. The emergence of the zero-gradient region indicates that the optimal control strategy switches between ON and OFF modes depending on the current state.

The OFF mode region appears if $\gamma$ is high and $N$ is moderately large (see the Fig. S1 in Sec. E of the Supplementary Material [27]), implying that the OFF mode is attributed to the difficulty of controlling an exponentially growing population. If both $\gamma$ and $N$ are large, the uncontrolled system has a stronger and less stochastic driving force towards $n_{A}=N$ as $n_{A}$ increases. Thus, the optimal control keeps the population away from $n_{A}=N$ by forcing it close to $0$, which is the ON mode around the peak $n_{A}^{*}$. Once a large fluctuation drives the population to leave the peak, the additional cost to bring it back to the peak does not pay for its success rate, leading to the OFF mode, i.e., do nothing, in the intermediate region. Finally, the control turns ON again near $n_{A}=N$ to hang on there. Such a spontaneous emergence of hierarchical control may not be identified within the LQG approximation, highlighting the importance of taking into account the unique properties of RN.

Lastly, we apply our framework to epidemic problems. We use the stochastic SIR model (Fig. 1 (c)) in which the population is divided into three classes: susceptible ($S$), infected ($I$), and recovered ($R$), i.e., $X=\{S,I,R\}$. The uncontrolled process is described by the following infection and recovery reactions:

and we assume mass-action-type kinetics, i.e., $h_{1}(n)=N^{-1}n_{S}n_{I}$ and $h_{2}(n)=n_{I}$. The total population size $N:=n_{S}(t)+n_{I}(t)+n_{R}(t)$ is a conserved quantity of the system. This model has absorbing states $\mathcal{N}_{\mathrm{abs}}=\{n\in\mathcal{N}|n_{I}=0\}$ and any stochastic solution $n(t)$ eventually reaches one of these states as $t\rightarrow\infty$. When the ratio $R_{0}:=k_{1}^{0}/k_{2}^{0}$ between reaction rate coefficients is high, the infection spreads rapidly in the population, and the state $n(T_{\mathrm{exit}})$ at the terminal time tends to have a small number $n_{S}(T_{\mathrm{exit}})$ of susceptible and a large number $n_{R}(T_{\mathrm{exit}})$ of recovered people. The number of recovered people at the end is equal to the total number of infections during the epidemic, which is known as the size of the epidemic [47].

The goal of control is then to minimize the size of the epidemic or, equivalently, to maximize the number of susceptible people at the end of the epidemic. Let us formulate it as the first exit problem where $U_{\mathrm{exit}}(n)=n_{S}$ and $U(n)\equiv 0$, i.e.,

We numerically calculated the value function and the optimal reaction rate using Eqs. (20) and (13) for $N=50$. The results are shown in Fig. 4. The value decreases as the number of susceptible people decreases and as the number of infected people increases (Fig. 4 (a)). The optimal control of infection and recovery rates plotted in Fig. 4 (b) and (c) indicates that strong control to reduce the infection rate and to boost the recovery rate is encouraged when the number of susceptibles is large and the number of infected is small. However, if many people are already infected, strong control is no longer encouraged. Instead, almost no control over the infection and recovery rates becomes optimal, as indicated by the white regions in Fig. 4 (b) and (c). This means that almost the entire population will eventually become infected in this situation, no matter how optimally the rates are controlled, and that further investment in the control cost is not worth the potential gains.

In particular, we can observe a sharp transition between strong control and no control when $R_{0}$ is high. Figures 4 (d) and (e) show the value function $V(n,1)$ and its gradient $\overline{\nabla}_{s_{2}}V(n,1)$ on the line $n_{S}=20$, which determines the optimal recovery rate $k^{\dagger}_{2}(n)=k^{0}_{2}\exp(\overline{\nabla}_{s_{2}}V(n,1))$. The value functions are approximately piecewise-linear functions $V(n,1)\approx\max\{0,a_{0}n_{S}+a_{1}n_{I}+a_{2}\}$, so the gradients and the optimal rates are approximately piecewise-constant. This transition leads to the mode switching of optimal control, which is similar to the case of the problem of maintaining diversity.