Reinforcement Learning for Control of Non-Markovian Cellular Population Dynamics

## Abstract

Abstract Many organisms and cell types, from bacteria to cancer cells, exhibit a remarkable ability to adapt to fluctuating environments. Additionally, cells can leverage memory of past environments to better survive previously-encountered stressors. From a control perspective, this adaptability poses significant challenges in driving cell populations toward extinction, and is thus an open question with great clinical significance. In this work, we focus on drug dosing in cell populations exhibiting phenotypic plasticity. For specific dynamical models switching between resistant and susceptible states, exact solutions are known. However, when the underlying system parameters are unknown, and for complex memory-based systems, obtaining the optimal solution is currently intractable. To address this challenge, we apply reinforcement learning (RL) to identify informed dosing strategies to control cell populations evolving under novel non-Markovian dynamics. We find that model-free deep RL is able to recover exact solutions and control cell populations even in the presence of long-range temporal dynamics.

## 1 Introduction

In order to survive, organisms must adapt to unpredictable environmental stressors occurring over diverse timescales. As a result, biological systems display remarkable adaptive capabilities, making it exceptionally challenging to control them through environmental modulation—an open and significant question in the physics of living systems. Two examples of adaptive systems with great clinical importance are cancer cell resistance to chemotherapy [Vasan2019], and bacterial resistance to antibiotics [Blair2015]. In both settings, constant application of a drug does not typically result in population extinction, as drug application also drives a certain fraction of the population to alter their phenotypic state to become drug-resistant. This new resistant state can persist even after drug removal and across cell lineages, thus encoding a memory of the stressful environment which allows populations to more quickly adapt if the drug is reapplied [Harmange2023, Mathis2017, Banerjee2021a].

Previous work [padmanabhan2017reinforcement, engelhardt2020dynamic, Gallagher2024, fischer2024minimising] has studied temporal drug dosing protocols to slow or prevent adaptation in various cancer or bacterial models using different methods. However, these models fail to account for the variety of adaptive timescales present in real biological systems. The presence of such memory effects greatly complicates the control problem, making the study of more realistic models imperative. Thus, to study the control of such memory-based adaptive systems, in this work we introduce a novel population model exhibiting phenotypic plasticity with non-Markovian dynamics. We couple insights from control theory with deep reinforcement learning to discover interpretable treatment protocols which successfully prevent proliferation.

## 2 Proposed Model and Approach

Figure: Figure 1: Left: Interaction loop between RL agent and environment (Sec. 2.2). Middle: depiction of the phenotypic switching model (Sec. 2.1). Right: Effect of learned policy on resistant fraction. for different memory strengths.
Refer to caption: x1.png

### 2.1 Non-Markovian Phenotypic Switching Model

To model the treatment response of an adaptive cell population, we use a general phenotypic switching model which captures the time evolution of a susceptible subpopulation, with size $S(t)$, and a resistant subpopulation, with size $R(t)$. Phenotypic switching models have been successful in describing a wide variety of biological scenarios, including development of persister cells and resistant cells in bacteria, and development of drug resistance in cancer cells [Balaban2004, fischer2024minimising, Witzany2023, kumar2019stochastic, Kratz2024]. In our model, susceptible cells with a net growth rate $\kappa_{S}(u)$ ($\kappa_{S}(u)<0$ corresponds to net cell death) switch to a resistant state at a rate $\alpha(u)$, where $u$ is the drug concentration. Similarly, when the drug is removed, resistant cells with a net growth rate $\kappa_{R}(u)$ switch back to the susceptible state at a concentration-dependent rate $\delta(u)$ (Fig. 1, middle). All growth and switching rates are a function of drug concentration normalized by the maximum allowable dose, thus $u\in[0,1]$. The time evolution of the size of each subpopulation, $\textbf{x}(t)=[S(t),R(t)]^{T}$, is then given by the dynamical system:

$$ $\dot{\textbf{x}}(t)=\textbf{f}(\textbf{x}(t),u(t))=\textbf{A}(u(t))\textbf{x}(t)$ (1) $$

where the state-transition matrix is given by:

$$ $\textbf{A}(u)=\left[\begin{array}[]{cc}\kappa_{S}(u)-\alpha(u)&\delta(u)\\ \alpha(u)&\kappa_{R}(u)-\delta(u)\end{array}\right]\;.$ (2) $$

We choose $\kappa_{S}(u)$ and $\kappa_{R}(u)$ to decrease and increase with drug concentration, respectively. This parametrization corresponds to “drug addiction” behavior in the resistant subpopulation, a robust phenomenon which has been observed not only in cell culture [SUDA2012292, Sun2014, Moriceau2015], but in animal models [DasThakur2013] and in vivo [Seifert2016138, Dooley2016]. Similarly, $\delta(u)$ decreases with $u$ while $\alpha(u)$ increases with $u$, as drug application drives the population to become more resistant, while reduction in drug concentration causes the system to recover susceptibility (see Appendix A for complete model details).

Recently, cell populations of many types, including human cancer cell lines, yeast, and bacteria, have been shown to maintain a memory of past environments which facilitates adaptation to previously-seen stressors over many timescales [Shaffer2020, Harmange2023, Larkin2024, Wolf2008, Mathis2017]. To better capture this memory dependence on treatment response, we introduce a memory kernel into the previously-described dynamics (Eq. (1)), making them non-local in time. Specifically, we choose a fractional differential equation (FDE) formulation as a phenomenological way to introduce multiple timescales of adaptation, one which has been used successfully to model memory effects in other biological [Lundstrom2008], ecological [Khalighi2022], and physical contexts [Bonfanti2020]. With this addition the dynamics now become:

$$ $\dot{\textbf{x}}(t)=\textit{{F}}(\textbf{x}(t),u(t))=\int_{0}^{t}\frac{(t-\tau )^{\mu-2}}{|\Gamma(\mu-1)|}\textbf{f}(\textbf{x}(\tau),u(\tau))d\tau\;,$ (3) $$

where $\Gamma(\cdot)$ denotes the Gamma function, $\textbf{f}(\cdot)$ is given by Eq. (1), and here we introduce the parameter $\mu\in(0,1]$ which controls memory strength. A value of $\mu=1$ corresponds to the memoryless case (first order derivative), whereas smaller values of $\mu$ correspond to an increased influence of past states on the current dynamics (lower order fractional derivative).

We seek to obtain a temporal drug concentration protocol $u(t)$ which minimizes the growth of a population and thus choose the final cost as $C:=\log N(T)/N(0)$ over the time interval $T$, where $N(t)=S(t)+R(t)$ represents the total population.
Thus, we aim to solve the following control problem:

$$ $\min_{u}C(\textbf{x}(0),\textbf{x}(T;u(\cdot)))\;\;\text{subject to}\;\;\dot{ \textbf{x}}(t)=\textit{{F}}(\textbf{x}(t),u(t)),\;\textbf{x}(0)=\textbf{x}_{0} ,\;u(t)\in[0,1]\;,$ (4) $$

where $\textbf{x}(T;u(\cdot))$ denotes the state of x at terminal time $T$ subject to control $u$ from $0\leq t\leq T$. Interestingly, in our minimal model (Eq. (3)), as long as $\kappa_{S}(u)$, $\kappa_{R}(u)$, $\alpha(u)$, and $\delta(u)$ are monotonic functions, we can show that the optimal control solution follows “bang-bang” control, regardless of model parameters and memory strength (see Appendix A for details). Thus, the continuous state-transition matrix of Eq. (2) can be simplified into a discrete one with binary control inputs, without altering the optimal solution (Appendix A). Despite this simplification, model control remains difficult, as the number of times and duration of drug application must be optimized. Constant application of the drug at the maximum dose ($u=1$) results in cell adaptation and proliferation (Fig. 2): a highly suboptimal solution to Problem (4) and a catastrophic result in the clinical context. Previous work [fischer2024minimising] has shown that in the memoryless case ($\mu=1$), the optimal solution for this type of model requires an initial drug application phase, followed by pulsing between treatment and pause phases at a regular interval dependent on the model parameters. However, as seen in Fig. 2, we find that the addition of memory ($\mu<1$) renders the control strategy of the memoryless case ineffective, as cells which have previously encountered treatment switch to the resistant state faster upon subsequent applications. Furthermore, in clinical or experimental setting, obtaining the values which parameterize Eq. (3) is usually not feasible, thus obtaining the appropriate switching frequency through direct computation via optimal control (OC) theory becomes impossible. As a result, we seek to learn the optimal policy directly through experience using reinforcement learning.

### 2.2 Reinforcement Learning

Reinforcement learning (RL) allows an agent to learn a drug protocol through experience (Fig. 1, left), despite the non-Markovian dynamics and without access to the underlying environment-specific model parameters. To formulate the RL problem, we define the relevant characteristics as follows: State: The state ($s_{t}$) of the agent’s environment is a list of the last $K=5$ estimates of the instantaneous growth rate, $c_{t}={\Delta^{-1}}\log N_{t}/N_{t-\Delta}$, where $\Delta$ is the simulation time between actions. Thus, a history of past observations is encoded in the state vector, a crucial design choice if the agent is to learn control without a recurrent hidden state.
Action: As motivated in 2.1, we choose a binary action space $u\in\{0,1\}$ representing whether the drug is applied or not, as this is suitable to recover optimal control.
Reward: As we seek to solve Problem (4), the reward is simply the negative growth rate, $r_{t}=-c_{t}$. Notice that with this choice of reward function, the sum of rewards across a trajectory simplifies as $R_{0:T}=\Delta^{-1}\log N_{0}/N_{T}$, ensuring the agent’s objective is aligned with a reduction in total cell population.
Dynamics: The dynamics of the total cell population $N_{t}$ are governed by the dynamical system described in 2.1, initialized to be fully susceptible ($\textbf{x}_{0}=[1000,0]$). After an action is executed, the simulation (a numerical solution (^1^11The numerical solution of fractional differential equations requires some care; cf. [garrappa2018numerical] for details.) of Eq. (3)) is computed with the action (dose) fixed for $\Delta=0.01$ hours to compute the next state ($s_{t+\Delta}$).

We use our own implementation of Double DQN [van2016deep] based on open-source code for DQN [stable-baselines3] to train the agent in this environment. After an action is taken, it is stored in an experience replay buffer for later use. We parameterize the value function $Q$ with a neural network (with parameters $\theta$), and train via SGD by sampling mini-batches uniformly at random from the buffer. The loss function is defined as the Bellman residual loss – the squared difference between left and right hand sides of the following equation:

$$ $Q^{*}(s_{t},u_{t};\mathbf{\theta})=r(s_{t},u_{t})+\gamma\max_{u^{\prime}}Q^{*} (s_{t+\Delta},u^{\prime};\mathbf{\theta})\;.$ (5) $$

As common in value-based algorithms, we use an additional target network and use exponentially annealed $\varepsilon$-greedy exploration. Further details on reinforcement learning can be found e.g. in [sutton&barto, rainbow]. We tune over several hyperparameters (whose values we list in the Appendix) to optimize performance. To encourage the agent to reduce the cell population while decreasing runtime, we terminate the episode if the number of cells ever exceeds its initial amount, forcing the agent to initially apply a continuous dose. All code to reproduce our experimental results can be found at [https://github.com/JacobHA/RL4Dosing](https://github.com/JacobHA/RL4Dosing).

## 3 Results

Figure: Figure 2: Performance comparison of constant drug application, solution for the memoryless case, resistant fraction-based pulsing technique, and policy learned by RL. For the fraction-based policy, an optimal lower and upper bound for resistant fractions are found through sweeping (Appendix A.1). The RL policy is capable of controlling the cell population better than any other scheme.
Refer to caption: extracted/5918128/figures/all_policies.png

We first test DQN in the memoryless case ($\mu=1$), for which an optimal controller is known [fischer2024minimising]. In this case, the optimal policy can be derived based only on two consecutive resistant fractions. However, given only the growth rates, the RL agent is able to reliably recover the optimal policy with only two frames and without any access to underlying model parameters. It does so by first applying the drug and then pulsing regularly, in agreement with the OC solution, as seen in Fig. 2 (rightmost plot).
With confirmation that RL can find the optimal dosing strategy in the memoryless case, we turn to the more difficult memory-based dynamics ($\mu<1$). We do not have a solution to Problem 4 in this regime, so we compare to two baselines: The memoryless protocol and a modified version in which switching times occur when the resistant fraction reaches a threshold value. We note that these baselines can have arbitrarily small switching times, which is practically infeasible. Remarkably, we find that using a small but bounded $\Delta$, our learned policy outperforms both of these baselines (Fig. 2). In addition, our experiments show that decreasing $\Delta$ beyond a certain threshold does not considerably increase performance (Fig. 4). Interestingly, the learned policy reveals that the agent initially maintains constant drug application before transitioning to a memory-dependent pulsing protocol. The agent increases the frequency of pulsing throughout the trajectory until saturation at the maximum rate ($\Delta^{-1}$) (Fig. 3). For cases with memory, this increase in dosing frequency can be understood as the result of faster cellular adaptation back to the resistant state after multiple drug encounters. Thus, to maintain a susceptible population, the policy’s pulse frequency must continuously be increased to compensate for the memory-based adaptability. We also find that the agent maintains a lower average resistant fraction at higher memory values (Fig. 1, right).

## 4 Discussion

In this work, we study the control of a highly non-Markovian model of adaptive cellular growth dynamics using deep reinforcement learning. Although we focus here on a specific parameterization most relevant to cancer, we expect this methodology to be applied successfully to other scenarios, including resistance development in bacteria. We find that deep RL is capable of recovering the known optimal policy for the memoryless case and can successfully find a policy for memory-based systems which prevents proliferation, all without access to the underlying model parameters. We utilize frame-stacking [nature-dqn] to encapsulate the history of the agent’s trajectory, but in future work we will test the use of recurrent policies to capture more nuanced long-term effects and perhaps further improve performance.
We also plan to extend our deterministic framework to the stochastic setting, as population heterogeneity is known to further complicate the control process.
This work demonstrates the benefits of combining OC and RL, showing how their frameworks can be effectively integrated. This connection can be further developed by using the memoryless OC solution to enhance RL training through reward shaping [ng1999policy, adamczyk2023utilizing] or pre-training [uchendu2023jump] techniques.

Experimentally measuring the growth rate of a pathogenic population is easier than determining its resistant fraction. Surprisingly, we find that RL can learn successful policies directly from the growth rate making it a promising method for use in clinical settings.

## 5 Acknowledgements

JA would like to acknowledge funding support from the NSF through Award No. PHY-2425180; the use of the supercomputing facilities managed by the Research Computing Department at UMass Boston; and fruitful discussions with Rahul V. Kulkarni. JCK acknowledges support from the Computational Biology Department and from Shiladitya Banerjee in the Department of Physics at Carnegie Mellon University.

## Appendix A Optimal control of the phenotypic switching model yields a bang-bang control solution.

We seek to find a temporal drug protocol $u(t)$ which solves the control problem

$$ $\min_{u}C(\textbf{x}(0),\textbf{x}(T;u(\cdot)))\;\;\text{subject to}\;\;\dot{ \textbf{x}}=\textit{{F}}(\textbf{x}(t),u(t)),\;\textbf{x}(0)=\textbf{x}_{0},\; u(t)\in[0,1].$ (6) $$

We choose the terminal cost to be the logarithm of the relative growth of the population throughout the course of treatment, $C(\textbf{x}(0),\textbf{x}(T;u(\cdot))):=\log N_{T}/N_{0}$, where $\textbf{x}(T;u(\cdot))$ denotes the state of x at terminal time $T$ subject to control $u$ from $0\leq t\leq T$, and where $N(t)=S(t)+R(t)$. Using Eq. (1) and the fact that $N(t)=S(t)+R(t)$, the model dynamics can be rewritten as a single fractional differential equation, namely:

$$ $D_{0}^{\mu}\phi(t)=f(\phi,u)=(\kappa_{S}(u)-\kappa_{R}(u))\phi^{2}+(\kappa_{R} (u)-\kappa_{S}(u)-\delta(u)-\alpha(u))\phi+\alpha(u)\;,$ (7) $$

where $D_{0}^{\mu}$ denotes the Caputo fractional derivative of order $\mu$ starting at $t=0$ [garrappa2018numerical], and here we drop the explicit time dependence for notational clarity. This can equivalently be written as the continuous delay differential equation (we use this form in the main text):

$$ $\dot{\phi}(t)=\int_{0}^{t}\frac{(t-\tau)^{\mu-2}}{|\Gamma(\mu-1)|}f(\phi(\tau) ,u(\tau))d\tau\;.$ (8) $$

To relate the net growth rate to drug concentration for each subpopulation, we assume that both rates take the general form:

$$ $\kappa_{S}(u)=\kappa_{S}^{\mathrm{max}}-(\kappa_{S}^{\mathrm{max}}-\kappa_{S}^ {\mathrm{min}})g(u)\;,\;\;\;\kappa_{R}(u)=\kappa_{R}^{\mathrm{min}}+(\kappa_{R }^{\mathrm{max}}-\kappa_{R}^{\mathrm{min}})g(u)\;,$ (9) $$

where $g(u)\in[0,1]$ is a monotonic dose-response function which relates drug dose to net growth rate, $\kappa_{S}^{\mathrm{max}}>0$ denotes the maximum growth rate of the susceptible subpopulation in the absence of drug application, and where $\kappa_{S}^{\mathrm{min}}<0$ is the maximum death rate of the susceptible subpopulation caused by application of the maximum drug dose ($u=1$). Importantly, $\kappa_{R}^{\mathrm{max}}>0$ denotes the maximum growth rate of the resistant subpopulation, which occurs in the presence of the maximum drug dose, and $\kappa_{R}^{\mathrm{min}}<0$ corresponds to the maximum death rate of the resistant subpopulation, which occurs when the drug is removed. This parametrization corresponds to “drug addiction” behavior, a phenomenon observed in several types of cancers [SUDA2012292, Sun2014, Moriceau2015, Seifert2016138, DasThakur2013, Dooley2016]. Similarly, the switching rates can then be defined as:

$$ $\alpha(u)=\alpha_{\mathrm{max}}g(u)\;\;\text{and}\;\;\delta(u)=\delta_{\mathrm {max}}(1-g(u))\;,$ (10) $$

where $\alpha_{\mathrm{max}},\delta_{\mathrm{max}}>0$ denote the phenotypic switching rates of cells switching from the susceptible state to the resistant state, and vice versa.

Given these definitions, the Hamiltonian associated with this control problem is then

$$ $H(\phi(t),u(t),\lambda(t))=\lambda(t)f(\phi(t),u(t))\;,$ (11) $$

where the trajectory of the Lagrangian multiplier $\lambda(t)$ is the solution to the costate equation [Gomoyunov2023]:

$$ $\lambda(t)=-\frac{\partial_{\phi}C(\textbf{x}(0),\textbf{x}(T;u(\cdot)))}{ \Gamma(\mu)(T-t)^{1-\mu}}+\frac{1}{\Gamma(\mu)}\int_{t}^{T}\frac{\partial_{ \phi}\lambda(\tau)f(\phi(\tau),u(\tau))}{(\tau-t)^{1-\mu}}\text{d}\tau\;.$ (12) $$

Applying Pontryagin’s minimum principle, we obtain the resulting inequality

$$ $H(\phi^{*}(t),u^{*}(t),\lambda^{*}(t))\leq H(\phi^{*}(t),u(t),\lambda^{*}(t))\;,$ (13) $$

which along with Eqs. (7) and (11) can be used to obtain the optimal control policy:

$$ $u^{*}(t)=\arg\min_{u}g(u(t))B(\phi^{*}(t),\lambda^{*}(t))\;,$ (14) $$

where $B(\phi,\lambda)=\lambda\alpha_{\mathrm{min}}-\lambda(\kappa_{S}^{\mathrm{max}}
-\kappa_{S}^{\mathrm{min}}+\kappa_{R}^{\mathrm{max}}-\kappa_{R}^{\mathrm{min}}
)\phi^{2}-\lambda(\kappa_{R}^{\mathrm{max}}-\kappa_{R}^{\mathrm{min}}+\kappa_{
S}^{\mathrm{max}}-\kappa_{S}^{\mathrm{min}}+\delta_{\mathrm{min}}-\alpha_{
\mathrm{min}})\phi$. As long as $g(u)$ is a monotonically increasing function of $u$, then the resulting control is said to be “bang-bang”, where $u(t)$ only takes extreme values. The switching times between maximum and minimum values of $u$ is determined by $B(\phi(t),\lambda(t))$, the switching function, yielding the optimal control solution:

$$ $\displaystyle u(t)$ $\displaystyle=0\;\;\text{if}\;\;B(\phi(t),\lambda(t))>0,$ $\displaystyle u(t)$ $\displaystyle=1\;\;\text{if}\;\;B(\phi(t),\lambda(t))<0,\;\;\text{and}$ $\displaystyle u(t)$ $\displaystyle\in[0,1]\;\;\text{otherwise}.$ $$

In principle, the optimal control trajectory can be obtained through numerical integration of the model dynamics (Eq. (7)) along with the corresponding costate equation (Eq. (12)). In the context of bang-bang control on non-Markovian systems however, using this approach can be difficult, as it requires careful choice of integration technique and update rule. Thus, we turn to reinforcement learning.

As shown above, the optimal control solution to Problem (6) follows bang-bang control, regardless of model parameters. This allows the continuous model of Eq. (1) to be simplified to a discrete model with binary controls without altering the optimal solution. This yields:

$$ $\dot{\textbf{x}}(t)=\textbf{f}(\textbf{x}(t),u(t))=\ \biggr{\{}\begin{array}[] {rcl}\textbf{T}\textbf{x}(t)&\mbox{for}&u=1\;\text{(Treatment Phase)}\\ \textbf{P}\textbf{x}(t)&\mbox{for}&u=0\;\text{(Pause Phase)}\end{array}\ ,$ (15) $$

where now there are two state-transition matrices, given by:

$$ $\textbf{T}=\left[\begin{array}[]{cc}\kappa_{S}^{\text{min}}-\alpha_{\text{max} }&0\\ \alpha_{\text{max}}&\kappa_{R}^{\text{min}}\end{array}\right],\;\hskip 20.0000 3pt\textbf{P}=\left[\begin{array}[]{cc}\kappa_{S}^{\text{max}}&\delta_{\text{ max}}\\ 0&\kappa_{R}^{\text{min}}-\delta_{\text{max}}\end{array}\right]\;.$ (16) $$

We use this formulation of the environment when training the reinforcement learning agent.

### A.1 Obtaining the optimal solution in the memoryless case

Previous work [fischer2024minimising] has shown that the optimal pulsing protocol for the memoryless case requires an initial drug application phase until the resistant fraction reaches some upper threshold $\phi_{h}$, followed by a pause phase in which the resistant fraction decreases to some lower bound $\phi_{l}$. This is followed by repeated cycles of drug treatment and pause phases where the switching time occurs when the resistant fraction reaches $\phi_{h}$ and $\phi_{l}$, respectively. To obtain the optimal values of $\phi_{h}$ and $\phi_{l}$ for our specific model parameterization, we swept over values of $\phi_{h}$ and $\phi_{l}$ between 0.1 and 0.9, with increments of 0.04, selecting the values which yielded the highest net death rate ($\phi_{l}=0.48$ and $\phi_{h}=0.52$).

## Appendix B Reinforcement Learning Details

We adapted DQN from Stable-Baselines3 [stable-baselines3] with a Double DQN action selection rule [van2016deep]. We train the RL agent for $3\times 10^{5}$ total environment steps. We limited the episodes to be of length $10^{4}$ steps (corresponding to 100 hours in simulation). An MLP of fixed size ($2$ hidden layers with $64$ dimensions each and ReLU activation) was used to parameterize the $Q$-function.

Regarding the environment and training, we list several key implementation choices: At the beginning of an episode, the state is zero-padded to always be of length $K=5$. We experimented with adding a penalty for allowing the number of cells to increase beyond the initial amount upon termination, but found this was not necessary for successful training. Borrowing terminology from the literature on Atari environments [nature-dqn], we stack $K$ frames (previous cost values) to form the RL agent’s state vector. Although initially we let $K$ be a $\mu$-dependent hyperparameter, we found a constant choice of $K=5$ to work well across the values of $\mu$ studied.

We find exponentially decaying the exploration parameter $\varepsilon$ to work better than the typical l inear annealing (with constant, positive final $\varepsilon$) scheduling. We conjecture that this is because non-greedy actions can be quite detrimental (causing the environment to terminate), and exponentially decaying $\varepsilon$ ensures some exploration continues to occur but with increasingly fewer random actions. To ensure the agent is not overly myopic (especially for such long episodes) we found a large discount factor of $\gamma=0.999$ (corresponding to an effective horizon of $H=(1-\gamma)^{-1}=10^{3}$) to be helpful. When training the agent, we wait until the completion of one rollout episode, and take as many gradient steps as environment steps have occurred.

### B.1 Hyperparameters

We find that sweeping over a range of hyperparameters (as shown in Table B.1) did not have a significant effect on performance, though for reproducibility we list the final hyperparameters used (for $\mu=0.7$) below in Table B.1.

**Table 1: Finetuned Hyperparameter Values for Double DQN**
| Hyperparameter | Finetuned Value |
| --- | --- |
| batch size | $32$ |
| buffer size | $100{,}000$ |
| exploration rate | $0.05$ |
| frames stacked | $5$ |
| gradient steps (UTD / RR) | $1$ |
| learning rate | $3.60\times 10^{-4}$ |
| target update interval | $1,000$ |
| discount factor | $0.999$ |
| learning starts | $10,000$ |

**Table 2: Hyperparameter Sweep Ranges**
| Hyperparameter | Sweep Values |
| --- | --- |
| batch size | $16$, $32$, $64$ |
| exploration rate | $0.01-0.2$ |
| learning rate | $10^{-5}-10^{-3}$ |
| target update interval | $1{,}000,5{,}000,10{,}000,30{,}000$ |
| target Polyak averaging | $0.95,0.99,0.995,1.0$ |

## Appendix C Further Experimental Results

Figure: Figure 3: The learned policy shows a resemblance to the optimal memoryless strategy, with an initial constant application phase followed by a pulsatile phase. However, in the case of memory-based dynamics, the frequency of pulsing must be increased over time as discussed in Sec. 3. Since the policy is eventually limited by the simulation time, the pulsing frequency becomes bottlenecked by our choice of time discretization $\Delta$ after $\approx 20$ hours. Despite this, the policy is still able to perform well with rapid pulsing.
Refer to caption: extracted/5918128/figures/learned_policy.png

Figure: Figure 4: We find that the simulation time can have a significant effect on RL performance. For each choice of $\Delta$, we run a random sweep of size $30$ over various hyperparameters, selecting the highest-performing run for each $\Delta$. Since smaller values of $\Delta$ require longer compute-times for simulations, there is a tradeoff between the amount of time (also, the inverse of max pulsing rate, which may be more relevant in clinical settings) and the best performance. We have chosen to use $\Delta=0.01$ throughout.
Refer to caption: extracted/5918128/figures/dt_data.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: 1. How many random seeds/trajectories did you use in each experiment? What is the faded area in Fig. 4?

2. Have you tested how robust is the resulted policy to parameter noise of the model?
soundness: 3
strengths: ### Presentation:
The paper is well written, and for me, it was easy to understand, in addition, it cites many relevant papers across all fields (RL, control theory, medicine).

### Methodology:
I think this paper bridges some of the gap between real-life problems to RL; while it is more applicative in terms of RL, it suggests how to do it right and efficient. The analysis in Sec. 2 enables RL to solve the problem (still need to show it). The experiments are extensive and show results for a few settings, and how it affects the model results.
summary: The paper considers the problem of controlling cell-population by adjusting drug dosage. In line with medical research, the authors propose a dynamics model with memory (non-Markov) for this setting, based on ODEs, and derive the theoretical optimal controller for drug administration, that is found to be a "bang bang" controller. Using this insight, the authors formulate this problem as an RL problem, then use RL with reduced action-space to solve it. The method is evaluated against a few naive baselines on a simulated environment, showing both superiority over the baselines, robustness to noise and changes of the memory of the system.
weaknesses: ### Presentation:
In the current form, the paper is written as two papers: One that formulates the problem and another that utilizes RL to solve the problem. For my opinion, it should have a more solid story, which would require some rephrasing. Considering the title, I suggest noting the difficulties of using RL for the raw problem, and then present your analysis for the optimal control which provides a solution which implies a much simpler problem for RL.

### Experiments:
I think that the naive baselines are not enough here, maybe add an RL variant that operates over the continuous action-space, showing the importance of your optimal control analysis?
Also, there are missing details for the experiments -- for how long the RL agent is trained? what

### Contribution:
My concern here, is that the method is motivated by real world, but eventually tested in simulation. Which shakes the ground of this paper, as it all revolves around the suggested dynamics model. I think another step towards bridging RL and real-life should be added. Since it is less likely to perform experiments that involve real-life environment and control, I would suggest showing empirically that the proposed dynamics model could indeed predict cell population -- using real-life data.

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: Do the authors have any intuition on how RL would perform in a continuous drug control scenario? 

If the constructed model might not fully describe the underlying biological dynamics (cf. p.4 l 214f.), causing OC to be inaccurate, how does RL solve this issue? 

How are the objectives of maximizing net death and maintaining the stability of the subpopulations connected? 

Given that learning an RL policy in a specific real-world application would not be feasible (similar to OC), how well does the resulting policy generalize regarding the parameters underlying the assumed model (beyond the memory size)? How can overfitting to the chosen training specification be prevented?
soundness: 4
strengths: The paper is well-motivated and regards a relevant real-world application. It is generally well-written, explained, and easy to follow. It provides a formally sound theoretical connection of the problem at hand to optimal control theory, extended to the application of reinforcement learning (including realistically obtainable observations) to tackle current shortcomings. Convincing empirical results show improved adaptability and generalization to varying memory levels of the underlying model.
summary: The paper proposes to model the non-markovian process of adaptive cellular growth dynamics in bacterial and cancerous populations to learn optimal drug dosing using reinforcement learning. It shows successful training of an improved policy for memory-based systems compared to optimal control, which is robustness to observation noise and memory strength perturbations.
weaknesses: While, in general, the problem is well-explained and elaborated, an example or illustration might have been useful. Also, the figure placement and referencing should be improved to further ease readability. The overall objective of reducing the population by minimizing growth to negative growth could be stated more clearly in the former sections (e.g., 2.1). In Fig. 2 (right), a comparison to the memory-less or optimal control approaches used as a baseline later on could have been insightful. 

Minor comments: 

- p.2 l.89 no middle part in Fig. 2 (referencing the wrong figure?)
- p.6 l.308 RL uses gradient ascent on the reward, not SGD.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 4
questions: - Given that mammalian cancer cells are targeted as an application scenario which generally have a cell division rate of about 24h, how realistic are the assumptions in the paper made about the delta parameter and the control frequencies applied by the pulsating policies?
- Can you give more biological motivation for the experiment with the switching memory strengths in Fig. 4 (right)?
- How important is the switch to FQF from DQN in section 5? How well would DQN do on the more challenging problems considered there? How important is the use of Noisy Nets for exploration over just using annealed epsilon-greedy? These points are not examined, which makes an assessment of the justification for the methodological choices difficult.
soundness: 4
strengths: - Very clear exposition, introduction to the problem and model development
- Contribution of a novel model useful for RL-based optimization of drug dosing
- Useful insights from control theory that allow simplification of the problem
- Optimal solution given for the memoryless case which serves as a sanity check for the RL-based method
- Convincing results of the RL-based solution compared to three baselines
- Clinically highly relevant problem overall
summary: The paper addresses the challenge of controlling and driving cell populations, such as cancer cells, toward extinction in environments where exhibiting adaptability and memory of past stressors. This is particularly difficult in systems with unknown parameters or those displaying non-Markovian dynamics such memory-based systems. The authors propose a new memory-based model for switching population dynamics, develop an optimal solution for the memoryless case, and apply model-free deep reinforcement learning (RL) to develop effective drug dosing strategies for cell populations with memory. The proposed approach is tested by comparing its performance against known exact solutions in controlled scenarios, demonstrating that deep RL can accurately recover these solutions. Additionally, the approach is evaluated in more complex environments where the memory strength of cell populations varies. The results demonstrate the method’s robustness in managing dynamic, memory-influenced systems.
weaknesses: - Limited methodological contribution: author's mainly apply existing RL methods to a newly developed simulation model (I see the model itself and the insights gained about it as main contribution)
- Not clear whether assumed time-scales are realistic

---

confidence: 3
contribution: 4
initial_rating: 8
presentation: 4
questions: - Why choose DQN and not another algorithm? SAC and others can also be trained in discrete action spaces. It would be nice to have a notion of how DQN compares to other algorithms (in supp. materials).

My main questions have to do with the feasibility of the proposed technique:
- How feasible is it to measure the cell growth periodically in real life settings? I guess one cannot grow cells very frequently and the amount of input data could limit the efficiency of the method proposed. Also, how feasible is it to provide drugs at the proposed rates, and would tighter constraints actually impact the policy much (could it be that a tight constraint on the frequency of the drug administration actually drives the policy to another local minimum?)
- I did not really understand if the plan here would be to train an algorithm offline on simulated data and test it in the real world (sim2real) or if the purpose of this method is to show on simulated data that RL can be used to solve the problem but this will require training the algorithm on real data. Or is the goal to get a standardized drug administration schedule that would fit all cases? 
- If the goal is sim2real, some evidence or a path forward to the applicability of this technique would be a nice add-on.
- If the goal is online training, one could maybe consider offline reinforcement learning to help solve this problem?
- How can one assess how good and representative the growth model is? Is there any evidence that sim2real would work?
- Alongside the previous points: is 300K steps for training realistic in real life?
soundness: 4
strengths: I've found the paper compelling and engaging. It is well written, goes through the necessary background effectively, has enough math to thoroughly back its claims and build the necessary background.
summary: This paper studies how non-markovian cellular population dynamics can be chemically controlled through the lens of RL.
While prior work focused on markovian assumptions or unrealistic assumptions, this paper tries to apply deep RL (DQN-based methods) to assess how effective this would be to design drug delivery patterns that optimally control the pathogen cellular dynamics.
The results of the RL-based method are encouraging: the authors show that in the markovian case, the algorithm is capable of finding the analytical solution, and in the non-markovian case a good handling of the cellular growth rate can be achieved.
weaknesses: - The paper mentions that a state encoding the memory is required, transforming the non-markovian dynamic into a Markovian dynamic. The usage of recurrent networks (or perhaps other models like transformers?) could be used to solve that problem without encoding for these states. I find this idea compelling. I wonder if using an RNN is enough material for a follow-up paper, even though it is crucially missing from this one. Hence, not having it within this publication could actually be detrimental for the publicity around this publication altogether.
- Only one family of RL algorithms is being used. Although the authors show that it is robust enough to discover the solution of the problem, it'd be nice to see how DQN compares to other techniques.

Meta-review: This paper presents an approach using RL to control non-Markovian cellular population dynamics, specifically focused on drug dosing strategies for cancer and bacterial cells. The key findings are: (1) development of a memory-based model for switching population dynamics, (2) derivation of optimal control solutions for the memoryless case showing "bang-bang" control is optimal, (3) demonstration that model-free deep RL can recover these optimal solutions and generalize to more complex memory-dependent scenarios, and (4) evidence of robustness to parameter variations and observation noise.

(b) positive

- good theoretical foundation combining optimal control theory with RL
- practical relevance for medical applications
- good empirical validation including comparison to analytical solutions
- good ablation studies and robustness analysis

(c) Weaknesses:

- limited validation beyond simulation environments (real-world feasibility questions around measurement frequencies and implementation)
- Could benefit from more algorithmic comparisons (addressed in rebuttal)
- The connection between theoretical and RL sections could be stronger

The authors managed to address the reviewers' concerns during the rebuttal. I vote for acceptance.

Decision: Accept (Spotlight)