Estimation Algorithm for a Hybrid PDE–ODE Model Inspired by Immunocompetent Cancer-on-Chip Experiment

: The present work is motivated by the development of a mathematical model mimicking the mechanisms observed in lab-on-chip experiments, made to reproduce on microﬂuidic chips the in vivo reality. Here we consider the Cancer-on-Chip experiment where tumor cells are treated with chemotherapy drug and secrete chemical signals in the environment attracting multiple immune cell species. The in silico model here proposed goes towards the construction of a “digital twin” of the experimental immune cells in the chip environment to better understand the complex mechanisms of immunosurveillance. To this aim, we develop a tumor-immune microﬂuidic hybrid PDE–ODE model to describe the concentration of chemicals in the Cancer-on-Chip environment and immune cells migration. The development of a trustable simulation algorithm, able to reproduce the immunocompetent dynamics observed in the chip, requires an efﬁcient tool for the calibration of the model parameters. In this respect, the present paper represents a ﬁrst methodological work to test the feasibility and the soundness of the calibration technique here proposed, based on a multidimensional spline interpolation technique for the time-varying velocity ﬁeld surfaces obtained from cell trajectories.


Introduction
Recruitment of immune cells to a tumor is a key parameter in cancer prognosis and response to therapy and the complex relationship between cellular, noncellular components and secreted chemotactic factors plays an essential role in directing the migration of both activating and suppressive immune cell types. In recent years with the development of the highly multidisciplinary Organ-on-Chip field (OOC) [1,2], microfluidic technologies are employed as valuable in vitro platform tools to build tumor microenvironments, with modular degree of complexity and to visualize and quantify immune infiltration in response to anticancer therapies.
The work that set the stage for Organs-on-Chip (OOC) was published in 2010 [3] and thereafter OOC constituted a dynamic field of research, and substantial effort has been devoted to creating realistic mimics of different organs in recent years [4][5][6][7].
The coupling with live-cell imaging may enable extraction of single-cell tracking profiles which can be processed with advanced mathematical tools. In this context, the present study is inspired by the modeling of the complex mechanisms behind the cell dynamics and interactions between immune (ICs) and treated tumor (TCs) cells in microfluidic chips. In this framework, several studies [1,7,8] were conducted on immunocompetent cancer-onchips to assess the effects of therapeutic drugs on TCs and on the possible reactions of the immune system. One of the first studies addresses the role of formyl peptide receptor 1/annexin a1 axis in anti-tumor response to anthracycline-based chemotherapy [7]. Timelapse recordings were performed in a microfluidic platform designed for oncoimmunology [1] research to assess physical and chemical contacts between malignant and immune populations. Some of the earliest applications of microfluidic cell culture technology focused on modeling specific steps in the cancer cascade, including tumor growth [9] and expansion [10], angiogenesis [11][12][13][14][15][16], progression from early to late stage lesions involving an epithelial-mesenchymal transition [17,18], tumor cell invasion [19] and metastasis [20]. Although Organs-on-Chip represent an increasing in importance field of research because it allows experimentalists to have major control on the objects of investigation, resulting in more accurate measurements, some aspects still remain difficult to understand. In particular, a quantification of the average chemical gradients present in the chip environment is difficult to be estimated.
Motivated by these laboratory experiments, we present an in silico mathematical model to describe ICs migration in the case of an efficient interaction with TCs treated with chemotherapeutic drug.
The main goal of this work is to build a bridge between experimental data coming from Cancer-on-Chip experiment given as particles and macroscopic mathematical models, such as the one proposed in reference [21], in order to gain further insights on the dynamics and short-range interactions between cells. Indeed, the present works goes towards the direction of the mean-field limit to unveil the statistical properties of the model. Recently, one of the authors of the present paper has derived rigorously in Wasserstein's type topologies the mean-field limit (and propagation of chaos) to the Vlasov-type equation, in the framework of generalizations of the kinetic model given by Cucker-Smale dynamical system, see reference [22]. Here we propose a simulation algorithm based on a hybrid macroscopic-microscopic chemotaxis model able to reproduce the main features of phenomena observed in microfluidic chip, such as migration of ICs and short-range interactions with cancer cells representing the sources of chemical gradients.
Existing methods used in the literature so far were dealing with the description of particle trajectories and are essentially represented by statistics on cell trajectories, as in the work by Agliari et al. [23]. Another possibility, already applied in other papers dealing with particle trajectories, see for instance reference [24], is to compute the distance between immune cells and tumor cells across time. However, we would like to stress that here our main concern is to show the feasibility of the proposed approach and to create a connection with macroscopic modeling of the same problem, in order to see what happens when we go to the limit (corresponding to the situation in which millions of cells move in the environment). Our final aim is indeed to be able in the future to simulate the behavior of immune cells in the tissue of an organ. In this framework, a novel strategy for the estimation of model parameters to perform the validation against real data and the calibration of the mathematical model is introduced. Such strategy involves from one hand the analysis of ICs trajectories coming from experimental data to describe realistically the average speed and the pathways of immune cells and the creation of a synthetic dataset for evaluating the effectiveness of the calibration technique for the estimation of model parameters. Moreover, in order to determine the velocity distribution of ICs, a stochastic component is obtained analyzing real trajectories of cells falling in the area under examination and then it is added to the deterministic velocity field.
The modeling of dynamics and interactions in microfluidic chip environment was already studied in [18,[25][26][27] with particle models and with cellular automata models in [28]. In [21] a fully macroscopic mathematical (PDE) model was considered, both for the chemical gradient which is seen as an average field and for the density of immune cells. With the PDE model we were able to describe long-range interactions in the chip environment and a first estimate of the chemical gradient driving ICs movement was obtained in [29].
Here, instead, we are interested in describing efficient short-range interactions between treated tumor cells and immune cells. As an example of this phenomenon, we considered the study and the experimental setting in [7], where tumor cells treated with chemotherapy drug release a chemical stimulus sensed by healthy immune cells, thus promoting their migration towards the tumor cells. For the construction of the model, we used a discrete in continuous approach where we coupled a reaction-diffusion partial differential equation (PDE) describing the evolution of the average substances released by the TCs in the tumor microenvironment with a particle model, where every IC in the system was considered to be a single entity and provided with specific properties, as previously done in [30] in a different experiment. Here the hybrid approach was revealed to be crucial in providing a proper description of the multiscale phenomenon that with a traditional macroscopic model would not be possible to fully decipher. In particular, our model describes the dynamics of each IC by means of ordinary differential equations (ODE) in such a way that every cell can be followed individually. Since the motion of ICs is driven by chemotaxis, the migratory activity is regulated by a chemotactic term which allows them to sense the gradient of the chemicals in the tumor neighborhood and of specific forces such as adhesion-repulsion which establish between cells.
The calibration of the mathematical model is supported by the development of ad-hoc parameter estimation procedure based on an interpolation technique for the approximation of the velocity field of ICs across time, i.e., multidimensional spline method previously introduced in [31] and described in Section 4.1. Such procedure, applied to synthetic cell trajectories dataset, provides accurate estimates for the values of unknown model parameters and demonstrates to be a promising approach for the validation of the proposed model against real data.

The Geometry of the Microfluidic Chip and the Related Computational Domain
The immune-oncology chip designed for the experiment consists of three main culture chambers for plating adherent TCs and floating ICs connected by a bridge of microcapillaries allowing chemical and physical contacts. In Figure 1 a picture of the two boxes is shown. With regards to the dimensions, the capillaries have, respectively, width and length of 12 µm and 500 µm, while the height is of 10 µm; however, since in the video footage the experiment is recorded at a fixed height, the third spatial dimension in our framework is neglected. The cross-sectional dimensions of culture chambers are 1 mm (width) × 100 µm (height). Further details about the chip design are illustrated in [7]. We also underline that here we are considering a 2D culture in the liquid with dying cancer cells adherent to the glass slide and mainly static, and immune cells floating.
In the experimental set-up, the timelapse observation area comprises left TCs culture chamber and central one with loaded ICs. Our goal consists of modeling the migration of ICs towards the TCs taking into account the different forces acting on the cells. To achieve such a result, we neglected the migration of ICs in the right chamber and through the microchannels, and we only focused on motility patterns of infiltrated ICs in the TCs left chamber.
Moreover, in order to better analyze the short-range dynamics and interactions between tumor and immune cells, we restricted our study on a subarea of the left chamber of length and height equal to 200 µm. This area corresponds to the subset [400, 600] × [200, 400] µm, and we defined it as Ω = [0, L x ] × [0, L y ]. In Figure 2 we present the setting of the experiment at time t = 24 h, with a focus on the area under analysis. The main reason for considering only a portion of the chip is motivated by: • the reduction of the computational cost in view of the optimization procedure for the parameter estimation; • the focus on short-range interactions.
Moreover, Ω contains four treated TCs, and this makes the area a good representative of cell dynamics, since in the laboratory experiment there are about 60 TCs in the entire chamber.

Original Contribution of the Present Paper
To describe such complex mechanisms with a multiscale nature, the model here proposed belongs to the category of hybrid models, involving the coupling between macroscopic and microscopic models.
As mentioned above, we assume the presence of forces between cells, which are established due to the chemoattractant and due to the nature of the cells itself. For this reason, in the current work, we refer to the discrete-continuous model (1)-(2) proposed in [30] to describe the morphogenesis of the posterior lateral line system in zebrafish. Such a model includes a classical reaction-diffusion partial differential equation (PDE) used to describe the evolution of the concentration of chemicals released by TCs in the tumor milieu, coupled with ordinary differential equations (ODEs) for the inter-cellular dynamics of ICs, which establishes both due to the presence of the chemoattractant and to the forces generated between the cells. In particular, the use of a particle model permits the highlighting of the role of single cells in the space, to analyze short-range interactions and to attribute specific characteristics to cells behavior that a macroscopic model would not allow. As mentioned in Section 1.1, for computational reasons and for the datatype at our disposal, cell tracks extracted from video recorded at a fixed level on the z-axis, we neglect the third dimension and consider only the 2D case. However, the third dimension can be easily taken into account by the model, but we do not expect important changes in the overall dynamics.
Here, chemical signals are described as an average field by means of a reactiondiffusion equation with a possible source or degradation term given in S term. Equation (1) was endowed with Robin boundary conditions to state the exchange of chemicals between the internal and the external environment. The gradient field f influences the evolution of cell positions, according to a second order equation of the form: where X i , i = 1, . . ., N tot,I , is the position vector of the i-th IC, N tot,I is the total number of ICs, X = (X 1 , . . ., X N tot,I ) contains all the positions, andẊ = (Ẋ 1 , . . .,Ẋ N tot,I ) the velocities. The vector Y = (Y 1 , . . ., Y N tot,T ), i = 1, . . ., N tot,T , stands for TCs positions, which are taken as constants in the problem, since dying TCs do not migrate as time evolves. The function F includes several effects: from the detection of the chemical signal f (chemotaxis) to mutual interactions between ICs (adhesion and repulsion) and ICs and TCs (adhesion and repulsion). All these effects take into account a non-local sensing radius. The term µẊ i represents damping due to cell adhesion to the substrate. We remark that the use of second order equations to describe the cell motion was first used in the seminal work [32] in the framework of Newton's equation of motion for particles. Indeed, even in this case the acceleration is not so high, this formulation describes the effects of the presence of an external force causing velocities and direction changes. Moreover, the friction term for cells immersed in a fluid is also taken into account. In order to have a more complete model able to represent the randomness characterizing cell motility, we also introduced a stochastic component in the velocity field in the spirit of Langevin model, see the review [33] and references therein.
In this context, in order to better reproduce the phenomenon under consideration such as the presence of chemical stimulus of treated cancer cells, we applied the following modification respect to the model in [30]: • we model the presence of tumor cells, and we also take into account the repulsion forces to avoid overlapping. Eventually, a slight overlay may occur between tumor and immune cell in the case of close interactions. However, it does not seem to occur in the video footage; • a chemical gradient concentrated around sources represented by cancer cells is considered and diffused in the environment; • Robin boundary conditions for the inflow of chemoattractant in the area under consideration are applied and adjusted to drive cell migration in a diagonal direction, as observed experimentally; • the alignment effect between cells is discarded since in this context is not present; • a chemotactic sensitivity term-i.e., receptor saturation-is added in the drift term of the equation of particles motion, with the effect that chemotaxis of cells is reduced in areas of high chemoattractant concentrations; • a stochastic component in the particle velocities is added to have more realistic cell trajectories in terms of randomness.
The equations are solved in a square domain Ω, introduced in Section 1.1, subset of the chamber where the experiment was performed and the approximation was carried out using a classical central difference scheme in space and the Crank-Nicolson scheme in time, although in Appendix A we prefer to present the general form of the scheme based on the θ-method.
The main novelty of our approach here is mainly represented by: • the modeling study on ICs behavior by considering different scenarios, see Section 3; • the construction of a calibration algorithm to find a trustable estimate of the most significant parameters of the model to assess the effects of chemical gradients on the ICs movement, see Section 4.
In particular, such calibration algorithm-based on the mathematical model-is built to test the possibility of validating the model on real data provided by experimentalists. The main goal of the present work is indeed the introduction of a robust procedure for the estimate of model parameters by means of the paths taken by ICs in the chip subarea under examination. It is worth noting that this is not an easy task, since to succeed in the model calibration we need to extract a common behavior from the time-varying trajectories of immune cells located at different points of the examined area. For this reason, in this first work based on this framework we propose a calibration strategy taking into account this variability and applying it to a synthetic dataset, in order to assess the effectiveness of our methodology.
We underline that the synthetic dataset of ICs pathways has been created reproducing qualitatively the real trajectories extracted from the video footage of the experiment, in order to have a "realistic" dataset. In more detail, a set of trajectories of cells sharing an average behavior has been produced by suitably tuning model parameters to have an upper and lower bound of cell speeds taken from the experimentally observed ones. With these synthetic trajectories, we have computed ICs velocity field at every observation time as a surface produced by the model itself, supposing to have a single population of immune cells thus showing an average behavior. Such velocity field has been then approximated using a multidimensional spline interpolation technique, described in Section 4.1 to have smoothing in time and space on the dataset to be compared.
Finally, this average field has been used as our target solution to be used in the calibration procedure, as described in detail in Section 4.2. The methodology for the error quantification is based on the minimization of a cost functional depending on the difference between the target velocity field of cell in the whole domain and the velocity field obtained for another choice of the model parameters during the optimization procedure performed with a local search method and we compared them at every time step. In addition to this, we also studied the introduction of further terms in functionals with the aim of improving the optimization results and show the soundness of our algorithm representing a starting point for investigations on experimental data in the next future.

Main Contents and Plan of the Paper
We introduce a hybrid model composed of a reaction diffusion equation, describing the time and space evolution of signaling substances released by treated TCs in a subarea of the main chamber, and of an agent-based model made up of second order differential equations for each IC. The boundary conditions associated with the partial differential equation are of Robin type, so that we can have a control on the amount of chemical exchanged with the surrounding microenvironment. In addition, based on experimental observation, an ad hoc parameter estimation technique was developed. To summarize the main contents of the present work, the mathematical issues faced in this study are identified into: • the development of a mathematical model describing the behavior of ICs in shortrange interactions in the microfluidic chip environment; • the development of ad-hoc parameter estimation techniques for time-varying velocity fields.
The plan of the paper is as follows. In Section 2 we describe the biological framework that inspired our study, and we introduce the mathematical formulation of biologically inspired models, and we present the adopted model. Section 3 is devoted to the modeling study on the scenarios suggested from the observation of the dynamics in the laboratory experiment performed with the numerical simulation algorithm reported in the Appendix. Section 4 contains the parameter estimation techniques here developed and the results obtained with the calibration algorithm. A sensitivity analysis is also performed and reported in the Appendix. Finally, in Section 5 we discuss the presented results, and the future aims of our work.

Biological Framework
The OOC technology allow the design and recreation of more sophisticated in vitro cellular microenvironments under physiological or pathological scenarios; as potential candidates for better prediction of human responses than animal testing cannot. This success is also due thanks to compatibility with a variety of microscopy techniques including live-cell high-content imaging and to the wide variety of cells and tissues that the chips can host, see [1].
One of the most challenging scenarios for the application of these devices is represented by cancer-immune relations due to very complex and not still completely discovered signaling modalities between ICs and TCs. Some attempts have been presented with the aim of modeling the effect of drugs on TCs, see for instance [34,35] where the microfluidic devices helped to control and evaluate the magnitude of cancer responses.
In particular, here we refer to the study in Vacchelli et al. [7], where timelapse imaging of microfluidic co-cultures were performed to investigate the motility patterns and crosstalk between ICs and TCs in the context of chemotherapy-induced anticancer immune responses.

Setting of the Laboratory Experiments
Details on the in vitro microfluidic experiments, type of cells involved and loading procedures can be found in [7,36]. Briefly, human breast cancer cells were previously treated with anthracycline-based chemotherapy and were cultivated in the left chamber. This treatment triggers the process of immunogenic cell death [37], according to which TCs release danger signals. The right chamber is loaded with unlabeled human peripheral blood mononuclear cells (PBMCs) from healthy donors (with normal expression of FPR1). PBMCs start migrating biased by the detection of chemical signal produced by TCs and after crossing microchannels enter in contact with dying cancer cells, engaging in stable interactions leading to TCs killing events.
In the chip, the chosen culture medium is neutral, which means that no exogenous substance is introduced. Timelapse imaging was performed using a microscope placed directly inside the CO 2 incubator for all the duration of the recordings. Images were taken every 2 min over a period of 72 h of migration. Immune cells tracks (<400) were extracted in left chamber (interval time: 24-48 h) using TrackMate plugin(ref) available in the Fiji/ImageJ software (https://imagej.nih.gov/ij/, accessed on 1 August 2020).
Our aim consists of designing a basic mathematical model able to capture the main aspects related to the migratory movement of the ICs with respect to TCs. Therefore, we only focus on the case of efficient interaction between dying cancer cells exposed to ICD inducers and immune cells from healthy donors. Concurrently, we want to develop an algorithmic calibration procedure to derive model parameter values as close to empirical observations as possible.

Mathematical Framework
In this context, it is important to remark that late years have experienced an increasing interest in the direction of developing techniques to combine experimental data and mathematical models, in order to produce systems, i.e., in silico models, whose solutions could be as close as possible to the experimental outcomes. Indeed, the success of informed models is mainly due to the consistent improvements in computational abilities of the machines and in imaging techniques that allowed a wider access to data.
The involvement of the immune system in all stages of the tumor life cycle, including prevention, maintenance and response to therapy is now recognized as central to under-standing cancer development from a systemic point of view. The increasing availability of experimental data and of treatment options, have represented a breeding ground to construct always more precise models. From a mathematical point of view, equations of different nature can be considered, depending on the type of analysis to be carried out. In macroscopic models, consisting of PDEs, any reference to the single constituents, the cells, are neglected, and macroscopic quantities such as the average cell densities are taken into consideration, see for instance the classical models [38,39] and their application to the immunocompetent microfluidic chip experiments in the recent work [21], allowing the simulation of long-range dynamics of immune cells driven by the chemical gradients secreted by cancer cells in the environment. The evolution of chemical stimuli in the tumor microenviroment can be also described by means of ODE models, see for instance the work [25] consisting of three ODEs that model the dynamics of effector and tumor cells and the cytokine IL-2, is one of the first to describe tumor-immune cells interactions. Another model in the same direction is due to Lee et al. [18], where changes in velocity and directionality of ICs are modeled with a system of ODE considering the combined effects of the interstitial flow, the tumor mass, the cytokines IL-8 and CCL-2, and the related receptors. Moreover, ODE models describing cell responses heterogeneity to death ligand (cytotoxic drugs) observed in laboratory experiments with single-cell techniques is proposed in [40]. In [28] the immune system is depicted as a network in the framework of cellular automata, where the nodes are molecules and cells and the arches connecting them are the influences each node has on others. Another kind of modeling consists of combining the macro approach expressed by a PDE, with the micro approach expressed by an ODE.
Hybrid models, developed including the macroscopic and microscopic scales, have been deeply studied in recent years, see [41][42][43][44]. In particular, in [41] a hybrid discretecontinuous model of metastatic cancer cell migration was developed, while in [45] a hybrid discrete-continuum approach was applied to model Turing pattern formation. In [43] a 3D agent-based model is proposed using an off-lattice approach to simulate vasculogenesis. For a comprehensive literature review of hybrid models sees the book chapter [46].
Here, in the framework of hybrid deterministic PDE-ODE models, we propose a discrete-continuous approach previously applied in [30], where immune cell motion is governed by ODE and the diffusion of chemoattractant released by cancer cells in the environment is described by PDE. The present work is complementary to the fully macroscopic approach proposed in [21], since it allows a zooming on immune system dynamics. Moreover, we enriched the equations of cell motion with a Brownian component to reproduce zigzag pathways of ICs.

The Model
Here we develop a model that can mimic cell migration under the effects of chemicals released by TCs on immune system dynamics. The considered TCs are experiencing apoptosis during which they liberate damaged-associated molecular patterns together with tumor antigens that elicit an immune response. For the sake of simplicity, the alarm concentrations are not differentiated, and we indicate them by means of the chemoattractant f , whose evolution is described by a PDE. At the cellular level, the model is discrete and includes the equation of the motion of each IC. We also underline that since we refer to experiment of treated cancer cells, according to the experimental setting we do not need to include TCs migration or duplication in our model. Let us summarize the main ingredients which compose the model.
For the IC motion we use a second order dynamic equation, which takes into account the forces acting on the cells that arise from the presence of the chemical signals and from the mechanical interactions between cells. The effect of cancer is described by a chemotactic term produced by the gradient of the f . The cell-cell mechanical interactions due to filopodia, consist of a radial attraction and repulsion depending on the relative positions of the cells, see [47] for experimental results in this direction. To describe this effect we refer to the mechanism introduced in [48]. The action of ICs regarding TCs is described by a radial repulsion term, depending on cells positions. Finally, we introduce a damping term, proportional to the velocities, which models cell adhesion to the substrate [49][50][51].
About the concentration of the chemoattractant, we associate a diffusion equation including a source term, given by the chemoattractant released by TCs, and a natural degradation term.
In the current context, the species under examination are TCs and ICs, but we underline that the setting can be made more complex with the introduction of a greater number of cell species and with the presence of an exogenous substance in the environment. Furthermore, in the present framework we are not taking into account the killing activity of ICs towards TCs as we already did in [21], since we do not have a quantitative data regarding the killing rate of immune cells in the chip environment. However, in the future also this feature will be included in the model. Let X i be the position of the i-th IC and f (x, t) the total chemoattractant concentration, the following equations are introduced: where γ, µ, ξ and η are given constants and F n (·), n = 1, 2, 3, 4 are suitable functions.
The term F 1 is the chemotactic term and therefore it is related to the detection of a chemical signal by i-th cell in its neighborhood and it is taken to be a weighted average over a ball of radiusR and centered in X i : where · being the Euclidean norm, is a truncated Gaussian weight function, and independently on i. A similar definition holds for the vector quantity F 1 .
In addition, the function F 1 contains the gradient of the chemical substance and a chemotaxis function χ( f ), representing the chemotactic sensitivity of ICs. Such chemotactic function, suggested in [52] by Lapidus and Schiller (1976) has the form: where k 1 represents the cellular drift velocity, while k 2 is the receptor dissociation constant, which indicates how many molecules are necessary to bind the receptors. Such function is known as receptor saturation function and has the effect of reducing chemotaxis of cells in areas of high chemoattractant concentrations. In the integral,R will be chosen larger than the IC radius R I , thus relation (5) describes a chemical signal that is sensed more in the center of the cell and less at the edge of the cell extensions.
Function F 2 Function F 2 includes a repulsion effect among TCs and ICs. In particular, repulsion occurs at a distance between the centers of two cells less than R 1 , where R 1 is defined as the sum of the two radii, R I and R T , of the immune and tumor cells, respectively. With this formulation, repulsion occurs when an IC and a TC start being effectively overlapped. We assume where the function K depends on the relative positions Y j − X i , namely: where ω rep TI is a constant. The quantity R 2 > R 1 indicates the radius of the ball centered in X i , in which the centers of the tumor cells can fall. Thus, we are considering all the tumor cells in proximity of the center of an immune cell.
Function F 3 Function F 3 includes adhesion-repulsion effects between ICs. In particular repulsion occurs at a distance between the centers of two ICs less than R 3 and takes into account the effects of a possible cell deformation. Conversely, adhesion occurs at a distance greater than R 3 and less than R 4 > R 3 , and it is due to a mechanical interaction between cells via filopodia. We assume where the function K depends on the relative positions X j − X i , namely: where ω rep , ω adh are constants. R 3 = 2R I will be chosen, so that the repulsion occurs when two cells start being effectively overlapped. Please note that function (11) gives a repulsion which goes as 1/r, r being the distance between the centers of two cells as we can find in [53,54]. The function (12) is the Hooke's law of elasticity. The last term in the first equation is due to the cell adhesion to the substrate (see for example [49][50][51]).

Friction
The addend µẊ i in the equation of the motion is due to cell adhesion to the substrate, with damping coefficient equals µ.
In the diffusion equation, only cancer cells are responsible for the production of the chemoattractant, so that where N tot,c is the total number of cancer cells, and In the previous formula R T is the radius of a cancer cell, considering that the source of chemoattractant is defined by the dimension of a single cell.

Initial Conditions
Initial data for Equation (4) are given by the position and velocity of each IC: Moreover, since TCs do not migrate and maintain their initial position Y during the whole time, the initial data for Equation (3) is provided by the chemoattractant produced by TCs at time t = 0:

Boundary Conditions
Now, let Ω = [0, L x ] × [0, L y ] our domain, for the chemoattractant we require the inhomogeneous Robin boundary condition: where b signals the similarity with the inhomogeneous Neumann boundary condition and regulates the exchange with the external environment.
The system of Equations (3) and (4) can be now rewritten as:

Stochastic Model
The model (16) is composed of deterministic equations, but a stochastic version for Equation (4) can be formulated. In fact, in recent years, several studies have shown that ICs exhibit an intermittent motion composed of a walking phase and of a zigzag phase. The walk is characterized by pause steps between the run steps, while during the zigzag, cells tend to turn away from their last turn directions and prefer to move forward in a zigzag manner [27,55].
This characteristic walk is here described by Brownian motion [56] as a first approach to the problem, revealing to be effective in reproducing the randomness of cell trajectories. The stochastic equation for ICs motion is: Equation (17) contains the stochastic contribution, where ψ(t) is a Gaussian white noise, and σ is the standard deviation of ICs trajectories. With this formulation, ICs are not only subjected to mechanical forces such as adhesion or repulsion, but also to random factors that might be related to unknown cell mechanisms. Some estimates of parameter σ based on experimental data are provided in Section 3.1.3.

Study on Different Scenarios: Numerical Tests
In this section, we look at some of the different scenarios the model can produce. We remark that the numerical simulations are performed in MATLAB c . The computational time for a simulation on the total number of frames T f = 681(N ∆ t = 680) takes about 260 s on an Intel(R) Core(TM) i7-9750H CPU @ 2.60 GHz 2.59 GHz.

Scenarios Representing Relevant Features of ICs Dynamics and Interactions
One of the motivations of the present study is to show the features of the proposed hybrid model in describing Cancer-on-Chip experiment. To this aim, we explore the different dynamics in three significant situations represented by the scenarios here proposed. During apoptosis TCs release alarm substances, which remain localized in proximity of the cells, that become attractors for ICs. In absence of an incoming chemical flow through the boundary, if the immune cells fall in the basin of attractions of the tumor cells, would remain trapped in their proximity all the time. From the experiments, however, we observed that ICs migrate also towards the bottom (where a reservoir of TCs was located) and to the left side of Ω, due to the presence of other cells in the surrounding microenvironment, generating a diagonal motion, as illustrated in Figure 2.
To allow ICs migration towards the sides of the domain, we added an incoming chemical flow through the boundaries, which caused changes in ICs orientation.
Moreover, to have control over the changes in orientation and direction of the migratory activity, it is necessary to find proper parameter values to attain a balance between the internal chemical concentration due to TCs and the concentration present at the boundary. For instance, if the chemotactic concentration at the boundary is higher than the internal concentration, ICs will sense the resulting gradient and move towards the boundary. On the opposite, ICs will not sense the gradient at the boundary. Thus, a balance between the internal and the boundary concentrations must be reached to have an IC attracted by the TCs and by the concentration on the sides.
For the numerical simulations we consider the domain Ω where four TCs are present (see Figure 2). The numerical positions of the TCs and the initial concentration of the chemoattractant released by the TCs are shown in Figure 3. We numbered the TCs from 1 to 4 to distinguish them. In our preliminary study, we also changed the number of TCs and their positions, with the aim of assessing the effect of their presence and their localization on the overall dynamics of ICs.
We assume the coefficients a and b in the Robin boundary condition (15) in such a way to have an incoming flux through the sides y = 0, x = L x and x = 0.
The laboratory experiment recording immune-cancer cells interaction has an overall duration of 24 h, but after data analysis we identified the presence of 97 ICs crossing the area in a time interval of 681 frames, which correspond to 1360 min (22 h 40 min) of video recording. To provide an idea of how ICs dynamics changes in time, we divided the time interval in six homogeneous parts, so that we can take a picture of ICs positions at six different instants T k for k = 0, . . ., 5 (T 0 = 1 frame, T 1 = 136 frames (4 h 32 min), T 2 = 272 frames (9 h 4 min), T 3 = 408 frames (13 h 36 min), T 4 = 544 frames (18 h 8 min), T 5 = 681 frames (22 h 40 min)).
The initial conditions for Equation (4) are assumed as the initial positions and speeds of the experimental cells.
In the following we consider three different scenarios that describe the interactions between TCs and ICs in a subarea Ω of the chip. Specifically, we consider:
Deterministic Motion including Cell Death; 3.
The first scenario is assumed to be a prototypal case, where the reciprocal mechanical forces between cells and the chemotactic stimuli are the only guiding forces of ICs migration, the second case takes also into account the possible effects of TCs death over ICs motion, and the third case contemplates the addition of the stochastic component Equation (A11), which modifies the deterministic IC trajectories.  In Figure 5 the time evolution of the migratory activity of ICs at six different times is depicted. At the initial time, only one IC enters the domain, while at time T 1 the number of ICs increases, and they are directed towards the tumor. As the time grows, most of the ICs approach the TCs and stay nearby, accumulating around them; while the others, guided by the inflow of chemical signal, move towards the left-bottom boundaries of the domain.

Scenario 2: Deterministic Motion including Cell Death
The biological experiment inspiring our work is related to the phase of immunogenic cell death, during which the TCs, previously treated with a drug, release alarm molecules sensed by the immune system. This latter reacts, attacking the tumor.
In the current scenario we focus on the description of the effects of cell death on the ICs dynamics. We remark that in this preliminary study, we did not insert a specific term in the model to describe the death of TCs due to ICs, but we directly turned off some TCs. Of course, this represents a simplification of what happens in reality, but is here applied for illustrative purposes. However, in the next future we will include in our modeling the death of cancer cells as a consequence of killing activity of ICs, similarly to the modeling proposed in [21]. Please note that in this simplified setting we let dead TCs physically remain in the domain, thus the effects of repulsion with ICs are still effective, but they do not release alarm molecules anymore. This provokes a decrease of the internal concentration due to the absence of previous sources. In the following, we turn off the third TC during the time interval [T 2 , T 3 ], and then we turn off also the fourth TC in the time interval [T 3 , T 4 ]. Figure 4 shows the level of the initial concentration (interval [T 0 , T 2 ]), while Figure 6 presents the evolution of the concentration after cell death. Comparing Figure 6 with Figure 4, we can observe how the overall concentration decreases after the two cells are killed. Figure 7 shows the ICs dynamics. At time T 2 both the third and the fourth cells are approached by ICs, while at time T 3 ICs start moving away from the third cell, which was previously killed and not releasing chemicals anymore. Between time T 3 and T 4 also the fourth cell is turned off and from time T 4 onwards it is not approached anymore. The effect of cell death of some tumor cells (the cells 3 and 4 depicted in Figure 3) on ICs dynamics is depicted in Figure 7. Please note that in the present case a different behavior respect to the one depicted in Figure 5 is observed, since here the accumulation around dead cancer cells does not occur, as expected.

Scenario 3: Stochastic Motion
Several chemical signals are produced at sites of dying TCs and then diffuse into the surrounding environment. ICs sense these chemoattractants and move in the direction where their concentration is greatest, therefore locating the source of the chemoattractants and their associated targets. If on one hand the deterministic model allows us to control the direction of motion of the ICs acting of the concentration of chemicals, on the other hand the complex migratory activity of ICs is better described by a stochastic model, that takes into account the variability of cells, over the preset deterministic mechanisms. For this reason, the features of this last scenario seem to be more suitable to qualitatively describe the real problem.
Estimates of standard deviations. Here we preliminary determine an appropriate value for the standard deviation σ in Equation (17). To this end, we consider the trajectory P i of the i-th IC, for i = 1, . . ., N tot,I and its corresponding smoothing S i , obtained with a moving average (the Matlab c smooth function), and we compute the variance of each IC trajectory: where T i is the number of frames the i-th IC spends on the domain Ω and we obtain the sequence Σ = {σ 2 i } i=1,...,N tot,I . Successively, we average the values of Σ and divide by 120 s, which corresponds to the time interval between two consecutive frames: obtaining a variance expressed in seconds. As a last step, we apply the square root to relation (19) and then obtain the standard deviation σ of real trajectories. From the analysis conducted on the total number of ICs, we observed cells with variances very far from the average, probably because of the presence of different cell species among the group of ICs. We remark that the immune cell population is heterogeneous, since it is composed by T-lymphocites, monocites, dendritic cells; however, a clear distinction among cell species is not possible currently. For this reason, in order to obtain homogeneous data, we neglect these cell trajectories in the computation of the standard deviation. Indeed, our mathematical model aims at describing an average behavior, and therefore the trajectories whose variance is too far from the average are discarded.
In the presence of more data, we will try to make a classification of cells based on their pathways (length, stops, tortuosity, etc.). This will be the subject of a future study.
The cleaned standard deviations computed from experimental trajectories are reported in Table 1.
ICs dynamics. Figure 8 captures ICs at six different times, while Figure 4 shows the evolution of the chemoattractant. Differences with the deterministic dynamics shown in Figure 5 can be highlighted. In the stochastic case, ICs are more spread.
The computations for this case are performed with a finer time spacing of ∆t = 10 s, with the aim of better describing by the model the complex mechanisms really happening to cells at smaller time scales, and then, the plots are taken every 2 min, which is the timeframe of video recordings.
As an example, in order to evaluate the effects of the stochastic component on the dynamics and establish a qualitative comparison with the deterministic model, we depict the trajectories of six randomly selected cells. In Figure 9 we assume the starting point of a given experimental cell and show the trajectories obtained with the ODE model (blue) and the SDE model (red), while the corresponding trajectory extracted from the video footage is plotted in black. As can be observed, the deterministic trajectories are smooth and tend, in most cases, to stop in correspondence of the TCs, showing that once the immune cells have fallen in the basin of attraction of the tumor, they tend to be trapped there. The stochastic trajectories, instead, have a more similar behavior to the real ones, showing that after some time spent near the tumor, the immune cells move to the boundary. So far the Brownian motion yields satisfactory results as shown in Figure 9, but we do not exclude in the next future to deal with more homogeneous and richer dataset to better identify the probability distribution and then apply a general Levý walk.

Parameters Estimation on Synthetic Data
In this Section we present the calibration algorithm for the estimation of some model parameters and the results obtained with it. All the model parameters are reported in Table 1. We underline that some of them are given from experimentalists since are taken from the laboratory experiment (size of the cropped area, number of TCs and ICs across time, standard deviation of cell pathways). Other parameters, such as for instance, the diffusion coefficient of chemoattractant, cellular drift velocity, receptor dissociation constant or the radii of cells are taken from the literature. For the model parameters whose values still need to be assigned, we assumed fixed most of them, such as the radii of action between cells or adhesion/repulsion coefficients among cells, setting them in a phenomenological way through the observation of their effects on the overall dynamics.
Then, in order to reduce the computational cost of the optimization procedure, we apply our calibration procedure only to 4 model parameters, i.e., the coefficient of chemotactic effect γ, and the chemical inflow, respectively, at the left, right and bottom boundaries b 1 , b 2 , b 4 , since we observed qualitatively a great effect of such parameters on the overall dynamics and, as a simplification of the model, we assume b 3 = 0. In particular, the coefficient γ affects the speed of ICs, while b 1 , b 2 , b 4 affect the directionality of immune cells.
We underline that the parameter estimation procedure here described is applied to the deterministic version of the model, i.e., (3) and (4). An extension of this procedure to calibrate the stochastic model (3)-(17) will be included in a future work.
In this methodological study we want to show the feasibility and applicability of such estimation technique for some significant model parameters of the hybrid model describing Cancer-on-Chip experiment. Our final goal is indeed not the exact fitting of model parameters, but is to show that we are able to qualitatively reproduce the dynamics observed experimentally, even if the quantification of the chemical gradients in the environment is currently impossible for biologists in this kind of experiments. In the next future, we will complement this methodology with macroscopic models deriving from them to be able to simulate immunocompetent behavior in organs and tissues, where millions of cells are present. To this aim, we applied the proposed procedure to a synthetic dataset produced by the model itself-but strongly inspired by experimental one-to assess the soundness of this strategy. In the future, we aim at having more data in such a way to apply the proposed strategy to a real biological dataset.

Data preparation.
To succeed in the calibration of model parameters, we need to extract a common behavior from the time-varying trajectories of immune cells located at different points of the observed area. Thus, here we propose a calibration algorithm taking into account this variability constructing a "realistic" synthetic dataset of ICs pathways that can reproduce qualitatively the real trajectories extracted from the video footage of the experiment. In particular, a set of trajectories of cells sharing an average behavior has been produced by suitably tuning model parameters to have an upper and lower bound of cell speeds taken from the experimentally observed ones. Then, ICs velocity field given by 2D surface is computed by the model itself at every observation time of 2 min (as the video footage timeframe). It is worth noting that the application of this strategy implies the assumption to have a single population of immune cells showing an average behavior.
The velocity field is then approximated using a multidimensional spline interpolation technique, described in Section 4.1 to have smoothing in time and space on the dataset to be compared.
Main steps of the calibration algorithm.
The calibration algorithm applied to estimate crucial model parameters can be summarized as follows: • the dataset representing the target solutions of our procedure computed from synthetic ICs trajectories is obtained by the model with fixed parameters; • a time-space approximation of the velocity fields is carried out with the spline technique presented in Section 4.1 and used as target solution of the calibration algorithm; • perturbed model parameters are used as initial guess of a global search algorithm minimizing the norm of the difference between the target velocity fields and the estimated ones, as explained in Section 4.2.
The results obtained with the procedure above are reported in Section 4.3.

Multidimensional Interpolation
The position and speed of the ICs, variable in time, provide punctual information about the velocity field in which the ICs are immersed. As a consequence, they can be considered to be training points for the calibration of an interpolation algorithm, able to provide information about the full space. Here we are considering time as the third dimension, being x and y the first two. As the interpolation scheme we are using a multidimensional spline, described in [31].
Each training point is supposed to influence the value of the interpolating function over a (limited) portion of the space in proximity of it. The degree of influence of the ith training point (w i (x, y, t)) is a function of time and space, driven by a compact-support function radial basis decaying at zero outside the influence area. Then, the interpolated value is obtained as a weighted sum of the values at the training points: where c(i) represents the value of the interpolating function at the ith training point and n is the total data-points available, corresponding to t 1 , . . . , t n discrete times. Specifically, we have two functionsṽ, i.e.,ṽ x (x, y, t) andṽ y (x, y, t). Moreover, the vector c(i) is different according to the velocities under exam, with training points corresponding alternatively to the xand y-velocities. Since different influence areas may overlap, the value of w i (x, y, t) needs to be adjusted to have a correct fit. To this aim, we need to solve a linear system where the n Equations (20) are collocated at the training points, generating an n × n system. The n coefficients h(i) are computed only once, and then applied in the interpolation procedure. In this specific case, we are using a linear function for the influence coefficients w i (x, y, t).

The Calibration Algorithm
With Formula (20) we have generated a velocity fieldṽ that from now on we indicate as V e to point out that the velocity field is produced with the punctual velocities v the immune cells assume at every time step, coming from experimental data (in this specific case from synthetic data). Later we have used this field to define an appropriate functional to be minimized.
To assess the goodness and soundness of our strategy, the analysis which follows is performed on synthetic data, i.e., numerical data produced by the PDE-ODEs system (3) and (4). Future work will be directed towards the application of this methodology to the experimental outcomes.
As a first step we have interpolated immune cell velocities v obtained from synthetic data to the end of generating a velocity field to use as a target, in correspondence of every frame. Thus, we have produced T f = 681 velocity fields V e x using the punctual velocities in x and T f = 681 velocity fields V e y using the velocities in y. In Figure 10 we present an example of the surfaces resulting from the interpolation at a certain time step. Figure 10A shows the interpolation of the x-velocities, while Figure 10B shows the interpolation of the y-velocities. The green points have coordinates (P x , P y , v x ) in (A) and (P x , P y , v y ) in (B). For more details about the interpolation technique see Section 4.1. Successively, we have launched the optimization algorithm and at every run we have computed the numerical solutions of the deterministic model with the approximation scheme described in Appendix A. We have then built a routine in which inserting the numerical positions we are able to find the corresponding velocities on the surface V e . The resulting interpolated velocities V n are then used to construct the objective function: In (21) we compare at every time step the punctual interpolated velocities with the punctual target velocities, i.e., the velocities used to generate the field. With θ we indicate the vector whose dimension corresponds to the parameter values we search.
In addition, we have also included a Tikhonov regularization term in the functional, as usually done for the regularization of linear inverse problems. In particular, we consider the following term to be added to the functional to be minimized: where θ 0 are the a priori estimate of target parameter values and θ are the parameter values we are optimizing. The constant λ is a regularization parameter that helps the algorithm in the search for optimal values of model parameters by reducing the number of local minima in the functional, see [62]. Tests for different values of λ showed better results of the calibration algorithm for λ 2 = 0.1. Thus, with relation (22) we want to reduce the error between the target values and the searched ones.
With functionals (21) and (22) we can define the minimization problem: where Θ is the space to explore to search the unknown parameter values. For completeness, we introduce another estimator we have used for our simulations, which is based on the idea of comparing the distances assumed by the tumor and immune cells at every time step. We indicate with C j the position of the j-th TC (the temporal indicator k is omitted since TC positions do not evolve in time). We call d k i,j the distance between the i-th IC and the j-th TC at time k: and then, fixed the i-th IC, we compare the distances d n,k i,j obtained from the numerical positions with the distances d e,k i,j obtained from the synthetic positions: , successively, we sum over all ICs: and then we sum over all times. To summarize, we have the functional: In conclusion, to evaluate the error committed by the optimization with respect to synthetic data, we indicate with θ 0 the target value and with θ * the corresponding optimized parameter and we define: the approximation error of each parameter. The minimization of the functionals above is performed with Particle Swarm Optimization (PSO) [63] as search method using the Matlab c toolbox. For the approximation of target and computed velocity fields resulting in a 2D varying in time surface we use a spline multidimensional interpolation described in Section 4.1.

Results on Parameters Estimation
Tests on the calibration algorithm are performed to assess the goodness and robustness of the methodology shown in Section 4.2. Please note that the total number of parameters to be assigned is 14, since the other parameters are given from the experimentalists or taken from the literature, as reported in Table 1. Of course, it is possible to make a parameter estimation of all of them. However, in the present work, in order to reduce the computational cost, we apply our calibration procedure only to 4 model parameters that we consider very significant since they strongly affect the dynamics of the ICs, i.e., γ, b 1 , b 2 and b 4 .
We recall the γ is the coefficient which enhances the effect of the chemotactic function (5), while b i , for i = 1, 2, 4 are the parameters that control the chemical inflow at the boundaries. In order to perform the tests, we produced synthetic solutions choosing specific values for the listed parameters: which correspond to those used to create the Scenario 1 (see Section 3.1.1). The synthetic solutions are used to generate the velocity fields V e to be used as targets.
To assess our strategy, we start testing one parameter and then adding one more at time. The range for the parameter variations is chosen by varying the initial guess by a percentage of ±50%. Specifically, we searched γ in I 50% = [3.5 × 10 −3 , 6.5 × 10 −3 ], b 1 in I 50% = [11,33], b 2 in I 50% = [6,18] and b 4 in I 50% = [9,27]. We also choose as initial guess a perturbation of the model parameters by +30%. We tested the algorithm with the following two different functionals: and Results with the functional J 1 are reported in Table 2, while results with J 2 are in Table 3. In each row of Tables 2 and 3 errors obtained with the parameter estimation algorithm are reported. Please note that errors are reported in increasing order respect to the number of parameters involved in the parameter estimation procedure, from the top to the bottom. We can notice that in both cases the errors in parameters estimation are quite low since their order of magnitude is around 10 −3 in the worst case, meaning the calibration was successful. We also underline that the norm of the total error of the functional is low, around 10 −7 , given by the J V component.   In conclusion, as a way to quantify the goodness of our reconstruction, we compute the percentage of cells passing through the cropped area of the chip representing the computational domain. In particular, we depict with histograms the statistics obtained from the trajectories reconstructed by the calibrated model and the statistics computed from the trajectories of the synthetic dataset, i.e., the target solutions of the calibration algorithm. The percentage of outgoing cells N k o over the number of cells N k p present in the domain Ω, is computed by the formula: for each timeframe k = 1, . . . , 680 with a time spacing of 2 min. Then, the optimized parameters (27) obtained with the calibration procedure described in 4.2 are used as the model parameters to be given in input to both the ODE model (3) and (4) and the SDE model (3)- (17) for the computation of the statistics on cell trajectories. More precisely, we perform this computation for the target case, i.e., the synthetic trajectories obtained by the PDE-ODEs system (yellow bars in Figure 11). Then, assuming randomly placed initial positions of ICs, we compute the same statistics for the PDE-ODEs system (red bars in Figure 11) and for the PDE-SDEs system (blue bars in Figure 11). In particular, the initial positions P k i , for i = 1, . . . , N tot,I are randomly perturbed in a range of (0, 5] µm. As can be seen in Figure 11, the statistics of target trajectories (ODE-target) and reconstructed trajectories (ODE-modified and SDE-modified) are quite similar in terms order of magnitude. Please note that the ODE-modified shows the ability to reproduce well also the timing of the exits from the domain, while a slight variability, as expected, can be observed in the SDE-modified case in terms of time occurrence of exits from the domain.

Conclusions and Future Work
In the present work we have developed a mathematical model to describe the interactions between ICs and treated TCs in microfluidic chip, together with a strategy to estimate unknown parameter values.
Regarding the modeling part, we have introduced a hybrid model, composed of a reaction-diffusion equation to describe the evolution of chemicals released by TCs, coupled with an equation for the motion of each IC, driven by the chemical gradient. The deterministic system represents a first microscopic model of ICs dynamics in a subarea of the chip.
First, we qualitatively analyzed the overall dynamics of ICs, varying important parameters of the system, such as the chemotactic coefficient γ and the boundary parameters. Indeed, according to their magnitude, the direction of the motion can be controlled and decided a priori. Then, in order to reproduce some interesting features characterizing cell behavior, we adjusted suitably the model parameters and we presented three different scenarios that can occur in the chip: • the deterministic scenario, with some ICs attracted by the tumor and others moving towards to the boundaries of the considered domain; • the cell death scenario, as a subcase of the previous one, obtained assuming two TCs have died after their interaction with ICs, causing changes in the internal concentration of chemicals thus affecting ICs dynamics; • the stochastic scenario, obtained adding to the equation of the motion a Brownian walk to mimic the randomness on ICs trajectories.
Regarding the stochastic scenario, it is important to highlight that the Brownian motion is not the most appropriate to describe ICs movement, but in the future we aim at substituting it by the Lévy walk. In recent works [64,65], it has indeed emerged that a heavy-tailed process is more efficient and realistic than Brownian motion. However, in the present work, the shortage of data at our disposal did not provide us information on the nature of the different ICs involved in the experiment. To this end, we are working on a classification strategy in a forthcoming paper to be able to identify the different categories of ICs and then differentiate the model parameters according to the corresponding cell type.
In addition to the pure modeling part, we have developed a model calibration procedure based on the comparison of the velocity fields related to ICs at every time step.
The calibration procedure with synthetic data revealed to be successful, thus representing a first step towards the model calibration on experimental data.
Finally, the main achievements of the present work are represented by: • the development of a simulation algorithm mimicking short-range dynamics of ICs in the neighborhood of TCs, as observed in Cancer-on-Chip experiment; • the introduction of an ad-hoc methodology for the calibration of model parameters based on the time-space approximation of synthetic velocity fields computed from immune cell trajectories.

Future perspectives.
Our future aim is to extend this framework to validate our model against the experimental velocity fields computed from the real IC trajectories extracted from the video footage of the experiment.
To face the problem of calibrating the model parameters against real data, we first must produce an interpolation model able to take into account the non-deterministic nature of the cell motions. A possible strategy is to produce a stochastic interpolator, able to determine the expected value and variance of the local speed: the final interpolated value will be then obtained adding a stochastic part, computed using the interpolated variance and the same probability density function of the experimental data, to the interpolated expected value. Using this approach, we have a non-deterministic interpolated value with the same statistical qualities of the interpolating dataset. Implicitly, we are considering as deterministic the local values of the expected value and of the variance of our experimental dataset, and we are also assuming to be able to compute the statistical distribution of the real data (otherwise, we can adopt a prescribed distribution. i.e., Gaussian). It is now evident that the comparison cannot be produced on a single path, but we need to observe the trajectories from a statistical standpoint, i.e., comparing the average path over many simulations, whose output is now non-deterministic. At the same time, in order to derive the statistical properties of the real IC trajectories and speed, a large number of experiments is needed, with a time resolution of the same order of magnitude of the time scale of the physical phenomenon under investigation. This implies a huge experimental and numerical effort, one or two order of magnitude larger than the present work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this section we discuss the numerical approximation scheme employed for Equations (3) coupled with (4) or (17)  , in this case if ∆t is the time step, t k will be the k-th temporal step, i.e., t k = k∆t, for k = 0, . . ., N ∆t . Note that the spatial steps were chosen as half radius of an IC: ∆x = ∆y = 2 µm, and the temporal step as the 1/6 of the video footage timeframe (of 2 min): ∆t = 10 s. With the notation u k n,m we denote the approximation of a function u(x, y, t) at the grid point (x n , y m , t k ).
Moreover, since experimentally it was observed that ICs leave the domain Ω, to manage the entrance and exit of cells and avoid numerical instabilities, we added a ghost grid to Ω d , where the cells lie after having left the main domain. To construct the grid, we considered the extended x interval [−L x * , L x + L x * ], with L x * > 0, with nodes x n = n∆x, for n = −N * , . . ., 0, . . ., N * , and the extended y interval [−L y * , L y + L y * ], with L y * > 0, with nodes y m = m∆y, for m = −M * , . . ., 0, . . ., M * . The equations were solved only on Ω d .

Appendix A.1. Discretization of the PDE
Regarding the approximation of the parabolic Equation (3). The right hand side is composed of the diffusion term, the source term, and the stiff degradation term −η f .
To eliminate this last quantity, we perform the classical exponential transformation: which leads to the diffusion equation with source for u(x, t): For this equation we apply a central difference scheme in space, i.e., the 5-point stencil for the Laplacian, and the parabolic Crank-Nicolson scheme in time.
The numerical scheme can be written as: where the second finite difference D 2 x u is defined as the central difference: and D 2 y u is defined analogously.

Appendix A.2. Boundary Conditions
Using the transformation (A1), the associated boundary conditions are: that we rewrite as ∂u ∂n + pu = q(t), with p = a D and q(t) = e ηt b D . Moreover, to distinguish the values of p and q on the different sides of Ω, we number them as follows: p 1 and q 1 are assumed on y = 0, p 2 and q 2 on x = L x , p 3 and q 3 on y = L y and p 4 and q 4 on x = 0. For the discretization of the boundary conditions, we use a central finite difference scheme. On y = 0 and y = L y , we have: for s = 2, 4 and k ≥ 0. The signs of r s and h s depend on the incoming/outgoing flow. For instance on y = 0, we have r 1 = −p 1 and h 1 = −q 1 , on y = L y , we have r 3 = p 3 and h 3 = q 3 , on x = L x , we have r 4 = −p 4 and h 4 = −q 4 and on x = 0 we have r 2 = p 2 and h 2 = q 2 in order to have an incoming chemical flux through the boundary. The previous discretizations (A5) and (A6) are used with the well-known ghost nodes technique [66].

Appendix A.3. Discretization of the ODE
The equation of the motion (4) is reduced to the first order system for i = 1, . . ., N tot,I . Equation (A7) is discretized with a one-step IMEX method, putting in implicit the term containingV i and in explicit the other addends, see [67]. Equation (A8) is solved with the forward Euler method. The two-dimensional integral in (A7) can be computed by a 2D quadrature formula, which due to the truncated Gaussian weight function w i (x) given in (6), is reduced to a sum of the discretized integrand functions on the grid points belonging to the ball B(X i ,R). For an integrand function g(x, t) holds: g(x, t)w i (x)dx ≈ ∇g(x n , y m , t k ) ≈ g k n+1,m − g k n,m ∆x , g k n,m+1 − g k n,m ∆y .
Equation (A7) is discretized as follows: and Equation (A8) is discretized as  (17) can be decoupled as follows: for i = 1, . . ., N tot,I . The discretization of Equation (A10) coincides with the one for Equation (A7), while Equation (A11) requires the application of the Euler-Maruyama method [68]. Equation (A11) can be written in the differential form and use dW(t) = ψ(t)dt where dW(t) denotes the differential form of the Brownian motion: where X i (t) is a one-dimensional Wiener process with drift V i and diffusion σ. This equation is discretized with the Euler-Maruyama scheme, which is the stochastic version of the deterministic Euler scheme. The increments of the Wiener process are defined as: The increment ∆W is a random variable with zero mean and variance equal to ∆t: ∆W ∼ N (0, ∆t), and with this increment we can construct approximations by drawing normally distributed numbers from a random generator. We approximate the process (A12) at the discrete time points t k , 0 ≤ k ≤ N ∆t − 1 by where ∆W = √ ∆tZ k , with Z k being standard normal variables with mean 0 and variance 1 for all k.