Towards open-ended evolutionary simulator for developing novel tumour drug delivery systems

Tumours behave as moving targets that can evade chemotherapeutic treatments by rapidly acquiring resistance via various mechanisms. In Balaz et al. (2021, Biosystems; 199:104290) we initiated the development of the agent-based open-ended evolutionary simulator of novel drug delivery systems (DDS). It is an agent-based simulator where evolvable agents can change their perception of the environment and thus adapt to tumour mutations. Here we mapped the parameters of evolvable agent properties to the realistic biochemical boundaries and test their efficacy by simulating their behaviour at the cell scale using the stochastic simulator, STEPS. We show that the shape of the parameter space evolved in our simulator is comparable to those obtained by the rational design.


Introduction
Precision oncology is a novel approach in tumour treatments that uses patient-specific information about the tumour to define the best possible targeted treatment (Ashely, 2016). Even when the right drug is identified for the right tumour profile, the problem remains of how to develop an appropriate drug delivery system that will pass all physiological barriers and reach the tumour in a high enough dose to be efficacious (Stillman et al., 2020). The sources of difficulty are numerous (Richardson and Caruso, 2020;Germain et al., 2020). Here, we will focus only on the problem of optimal parameter identification. Our question is how to identify an optimal set of parameters, such as binding affinity, uptake, and diffusion rate in order to reach desired treatment efficacy, having in mind that tumour can rapidly develop drug resistance. With the recent increase in clinical research on developing targeted therapies with multiple drugs to deal with the adaptable tumour, two modelling problems came into focus. The first one is how to design optimal treatment that includes several drugs (Tsompanas et al., 2020. The second, more general problem, is how to design treatment for the adaptable tumour even when future phenotypes are unknown. Recently we have developed a simple agent-based, open-ended evolutionary engine to automatically design combinatorial oncological treatments . There, we demon-strated that open-ended evolution coupled with multidimensional optimization can lead to the emergence of successful virtual therapies. Here we introduce several improvements described below.

Model description and simulation results after mapping parameters to realistic values
Details of the original model and the range of all parameters are given in Balaz et al. (2021) and here we only give a brief overview. In the model we previously introduced the simulated world is represented as a 2D grid where evolvable nano-agents (NA) learn to eliminate tumour cells (CC) while not harming healthy cells (HC). NA can move, have internal memory, and can learn about their environment. They can only "see" the contents of the grid position they occupy. When the NA "steps-on" a cell, it will check whether it already has it in its memory (by looking at the visible properties of the cells), and memorize it if not. If the memory is full -the new cell information will take place of the oldest memorized cell. At the beginning of the simulation all agent memories are empty. Interaction of NA agents with cell-agents is modelled as a Michaelis-Menten reaction network: where NA f is a free NA, R is a receptor at the cell surface, C is a NA-cell complex and NA i is an internalized NA. Probabilities of association (p a ), disassociation (p d ) and internalization (p i ) govern the dynamics of the NA-cell interaction. If NA encounters a "familiar cell", the probability of association with that agent is p a . Otherwise, the probability of association is reduced by multiplying p a with the "curiosity" factor. The local fitness value is assigned to each NA agent as a sum of killed CC minus the sum of killed HC. After each 10 time steps selection/mutation routine takes place. Mutable parameters are movement speed, p a , p d , p i and probability of killing (p k ). CCs can evade NA by mutating their visible properties by which NA recognize them.
To mimic tumour resistance, 10% of the total number of tumour cells were chosen randomly and assigned randomly- Cell 06 Cell 01 Cell 02 Cell 03 Cell 04 Cell 05 Cell 06 Free receptors NP-receptor complexes Internalized NPs Number of NA chosen resistance modifiers (one per cell). Resistance modifiers can change p a , p d , p i , and p k of NA that interact with them. Resistance strength is randomly chosen in the interval of 30 − 80% and is inherited after cell division.
The novel feature we introduce here is the split of the model into two modes: the learning mode and the simulation mode. In the learning mode, the goal is to find the best collective strategy to fight the tumour, regardless of the duration of the treatment. After interacting with a cell, NA agents do not perish but continue with movement. We ran the learning mode for 10, 000 time steps and the efficacy of top performing NAs is evaluated in the "simulation mode" and in the STEPS simulator. In the simulation mode we include the realistic temporal dimension based on the following assumptions: (i) one grid cell equals one cell agent, (ii) average tumour cell diameter is 10 µm, (iii) the diffusion coefficient range for nanoparticles in a fluid is in the range of 10 −10 cm 2 /s (Katayama et al., 2009;Pinheiro et al., 2007). Given that we are in 2D space, NA then needs t = (10 −3 ) 2 ÷ 2 × 10 −10 = 0.5 × 10 4 seconds to move from one empty grid cell to a neighbouring cell. As we fixed NA speed in an empty cell to 1, one time-step in the simulation mode equals 5000 seconds. To calculate injection dynamics of NAs in simulation mode, we took pharmacokinetics values for a typical anticancer drug Doxorubicin (Dawidczyk et al., 2017). There, a drug reaches maximum concentration in the tumour in 20 hours (∼ 14 simulated time steps) after the injection. For the following 80 hours (∼ 72 time steps) a drug concentration slowly declines until it reaches zero. Therefore, during the first 14 simulation steps, at each step we inject 1/14 of the predefined total dose. As expected, most efficacious option is when the maximum tolerated dose of NA is injected. This eliminated 7% of tumour cells after the single dose treatment (data not shown). To test how the spatial distribution of cells affects nanoparticle transport, we used the Stochastic Engine for Pathway Simulation (STEPS) (Wils and De Schutter, 2009) where we incorporates spatiality through the discretisation of the domain into tetrahedral well-mixed subregions. To normalize evolved NA parameters to realistic values we take the range of association rate constant (ka) to be between 10 4 and 10 6 [1/Ms] (Hauert et al., 2013). Since the time step is 5, 000 seconds, ka is 5x10 7 to 5x10 9 [1/M time step]. If we assume that one NA represents 10 5 particles (1.66x10 − 7 [M]), then the ka range is 8.3 to 8.3x10 2 [particles / time step]. The normalization factor is then 1.66x10 − 7 [M] / 5, 000 [s] = 1.2x10 3 . Therefore, for example p a = 0.3 translates to ka=3.61x10 2 [1/Ms]. Using the same approach we calculated the normalization factor for disassociation constant kd and internalization constant (ki) to be 2x10 −4 . In the simulations, maximum penetration depth is 6 cells Fig.1 and NA were able to kill all cells into which they were internalized.

Conclusions
At this stage it is hard to directly compare our results to clinically relevant ones due to the differences in tumour size scaling, drug-action mechanisms and influence of tumour heterogeneity on drug diffusion. Nevertheless, it is important to note that the distribution of NA parameters we obtained as a result of in silico evolution is similar to those used to synthesize nanoparticle-based DDS (which is result of highly skilled rational design) as well as to parameters obtained via deterministic modelling (Hauert et al., 2013). Therefore we believe that our evolutionary approach can generate useful leads in designing novel drug-delivery systems, especially since our engine is designed to deal with tumour changeability which is one of the main problems in long term efficacy of oncology treatments. Therefore, our ongoing work is focused on further closing the gap between our evolvable simulator and pre-clinical investigations, mostly by devel-oping closer mapping between the parameter space of the simulator and the corresponding biochemical counterparts.