Modeling Living Cells Within Microfluidic Systems Using Cellular Automata Models

Several computational models, both continuum and discrete, allow for the simulation of collective cell behaviors in connection with challenges linked to disease modeling and understanding. Normally, discrete cell modelling employs quasi-infinite or boundary-less 2D lattices, hence modeling collective cell behaviors in Petri dish-like environments. The advent of lab- and organ-on-a-chip devices proves that the information obtained from 2D cell cultures, upon Petri dishes, differs importantly from the results obtained in more biomimetic micro-fluidic environments, made of interconnected chambers and channels. However, discrete cell modelling within lab- and organ-on-a-chip devices, to our knowledge, is not yet found in the literature, although it may prove useful for designing and optimizing these types of systems. Consequently, in this study we focus on the establishment of a direct connection between the computer-aided designs (CAD) of microfluidic systems, especially labs- and organs-on-chips (and their multi-chamber and multi-channel structures), and the lattices for discrete cell modeling approaches aimed at the simulation of collective cell interactions, whose boundaries are defined directly from the CAD models. We illustrate the proposal using a quite straightforward cellular automata model, apply it to simulating cells with different growth rates, within a selected set of microsystem designs, and validate it by tuning the growth rates with the support of cell culture experiments and by checking the results with a real microfluidic system.

cellular Potts, Glazier-Graner, agent based, among others) 11,12 . These discrete models have some drawbacks when compared to continuum approaches, including computational cost for larger cell numbers and precise lattices and need for calibration upon macroscopic measurements. However, discrete models can be more easily fine-tuned by means of averaged measurements from controlled experiments, when the model parameters from continuum models are related to difficult-to-measure cell scale phenomena 12 . In this study we focus on modelling collective cell behavior by using discrete cell models, whose origins and applications to modelling cell colonies are detailed below.
Going to the origins of modern discrete modelling, cellular automata were developed on the basis of work by pioneers, such as Stanislaw Ulam and John von Neumann, as a collection of elements or cells defined upon a grid that evolves through time steps following a set of rules applied iteratively. Along the time steps, the state (i.e. colour or value, typically "0" or "1") of the cells within the grid changes according to the rules and to the previous states of neighbor cells 13 . Since the beginning, these models were conceived as possible simulators for biological systems and well-known examples of application appeared, such as Conway's game of life 14 , in which the cells upon a two-dimensional grid have two possible states, dead or alive, and in which cells survive, reproduce, die by under-or over-population, depending on the 8 neighboring cells or the previous generation. Apart from the initial game-like demonstrations, further studies led to verifying that extremely complex systems could be modeled by using cellular automata 15 .
More recently, in the specific area of modelling cell behavior, cellular automata have been used for modelling cell adhesion and proliferation; 16 for modelling migration, proliferation and differentiation 17,18 ; or, in connection with lattice-Boltzmann methods, to model multi-scale avascular tumor growth coupled with nutrient diffusion and immune competition 19 . As for other discrete cell models working upon lattices, the cellular Potts model 20 complements the lattice with an energy function or Hamiltonian that can be defined to control different cell behaviors, including migration, clustering and growth, and to add volume and surface constraints to the model. This approach has led to the implementation of CompuCell3D 21 , one of the most used software worldwide for modelling cells and their collective behavior, which has been employed for modelling cancer growth and invasion 22 , to simulate epithelial-mesenchymal transitions 23 , and also as educational tool for biomedical engineering degrees 24 , to cite just some examples selected among dozens of publications available in the CompuCell3D website (http://www.compucell3d.org).
In any case, all these discrete cell models used for predicting collective behaviors normally operate on infinite or boundary-less 2D lattices, hence modelling cell growth, migration, death and interactions in "Petri dish" like environments, with the same limitations as the use of Petri dishes for understanding in vivo performance using an in vitro approach. To the best of our knowledge, these cellular automata and cellular automata-like models have not yet been applied to modelling cells within lab-or organ-on-a-chip devices, which could support the systematic design and optimisation of these sorts of innovative biomedical devices and their steady regulatory compliance verification, by means of selected experiments and exhaustive simulations.
Consequently, in this study we focus on the establishment of a direct connection between the computer-aided designs (CAD) of microfluidic systems, especially labs-and organs-on-chips (and their channel structures), and the lattices for discrete cell modelling approaches aimed at the simulation of collective cell interactions, whose boundaries are defined directly from the CAD models. We illustrate the proposal using a quite straightforward cellular automata model, apply it to simulating cells with different growth rates, within a selected set of microsystem designs, and validate it by tuning the growth rates with the support of cell culture experiments and by checking the results with a real, although simple, microfluidic system. The materials and methods employed are detailed further on, before presenting and discussing the key results.

Materials and Methods
Creating the grid and boundaries for a cellular automata model from the CAD file. Throughout the study we use NX-8.5 (Siemens PLM Solutions) for computer-aided design purposes -mainly for designing the microfluidic systems and organ-on-a-chip devices-and Matlab (The Mathworks Inc.) for developing the code of the cellular automata model working upon such microsystems and performing the collective cell simulations. As mentioned earlier, a key objective of our model is to directly link the computer-aided designs of organ-on-achip devices with the lattices used for the cellular automata model. Accordingly, the channels and chambers of the microsystems should define the allowed cell positions limited by the vertical walls, which prevent cells from escaping the microsystems, as they operate typically closed by a microscope glass cover slip in real applications.
Such a connection can be performed in a quite straightforward fashion: Starting with the 3D CAD file of the microfluidic system, the model is positioned to show a top view, whose surfaces are painted using "paint operation" tools, so as to display the channels in white and the boundaries in black. A final conversion into. jpeg format enables direct import using Matlab as working lattice. The process has direct connections with the generation of digital masks, directly from CAD files, for the manufacture of microfluidic devices by mask-less UV-photolithography 25 . Figure 1 provides a simple example by showing the computer-aided design of a device, made of a single channel connecting an inlet and an outlet, and the obtained boundaries for the lattice-based simulation, in which the white pixels represent allowed regions for the cells (one cell-one pixel), while the black zone is prohibited and not considered for the simulation. As can be seen, each pixel or point of the lattice is defined by its position [X, Y] and by a colour in [R G B] format. The colour helps to distinguish between prohibited lattice points (black), empty lattice points (white) or lattice points with different cell types (i.e. green and blue), with dead cells (red) or cells affected by a drug (pink), to cite some options.
In order to avoid scale problems, the size of the mask should be carefully selected according to the real size of the microsystem and to the actual size of the cells under study, which defines each pixel, in our case corresponding to a square of some 20 × 20 μm 2 . Having a pixel per cell may be in some cases computationally expensive, so an alternative option would be to generate coarser lattices and to work with cell aggregates. Besides, in some cases a fine-tuning of the generated lattice, which derives from the digital mask created from the CAD file, is needed before starting the simulations. The reason is that sometimes the boundaries between the allowed (white) and prohibited (black) zones are not so sharp, due to exporting from CAD to a.jpeg file, and some slightly grey borders appear, which should be transformed to real white with y 0 = 1; d) models interactions between different cell types, defined with different colors, enabling the simulation of invasive cells (i.e. blue) that convert surrounding healthy cells (i.e. green) in invasive cells, as in Fig. 2b; and e) models the addition of a substance, drug or reactive, which diffuses following a proliferation equation similar to the one already discussed but with another dynamic, as it is defined to interact with a specific type of cell line. Cell migration is not considered in this first model implementation, as the cells we would be testing (see Section 2.4) are of www.nature.com/scientificreports www.nature.com/scientificreports/ adherent nature, but could be of direct implementation to achieve a more universal simulator, which could be also complemented with energetic functions similar to those used in the cellular Potts model and related ones. For the purpose of connecting CAD models with cellular automata with defined boundaries for simulating collective cell behaviors within organ-on-chip devices this preliminary approximation may prove sufficient. A simulation starts by defining the initial positions of the microsystem where the cells are placed, either by modifying pixel color of the desired positions by writing within the lattice matrix or by clicking in the ad hoc developed interface to position cells of different types and eventual substances, reactives or drugs. Each step goes on by evaluating the allowed positions in loops along x and y directions, by leaving white positions without neighboring cells white, by evaluating the probability of death of each cell (and eventually transforming those dead into red), by making healthy cells proliferate according to scheme 2a and by taking account cell-cell interactions, according to scheme 2b.

Set of computer-aided designs and prototypes of microsystems for testing the model. A set
of computer-aided designs was used to test the proposed simulation process. The models are shown in Fig. 3 and include: the already mentioned simple channel design (3a); a design with a central chamber and multiple radial channels, adapted from previous designs by our team and conceived for studying cell communication (3b) 26 ; a multi-chamber organ-on-chip with a central vascular channel, aimed at studying metastasis (3c); and an organ-on-chip device designed for studying interactions of the blood-brain barrier also adapted from previous designs by our team (3d) 27 . In our opinion, the selected set of designs provides enough versatility to test the cellular automata in connection with real examples of labs-and organs-on-chips.
Cell culture experiments for adjusting and validating the simulation. Two different cell lines N2A (ECACC 89121404) and MC3T3 (ECACC 99072810) were used, in the facilities of the UPM Centro de Tecnología Biomédica, which helped to adjust the growth rate of different cell types by means of Petri dish cultures and to incorporate these growth rates into the cellular automata model, running upon the multi-chamber microsystem, in which finally the cells were also cultured to validate the simulation.
Both types are adherent cells, in connection with our cellular automata model, in which the cells live, die or proliferate, but do not migrate. The N2A cells are a mouse neuroblastoma cell line with a neural and amoeboid stem cell morphology, which can differentiate into cells with features of neurons 28 , while the MC3T3 is an www.nature.com/scientificreports www.nature.com/scientificreports/ osteoblast precursor cell line derived from mouse calvaria 29 . Both cell types are cultured using the Dulbecco's Modified Eagle Medium (DMEM) enriched with 10% fetal bovine serum (FBS), 2 mM glutamine 2 mM, and 1% penilicin-streptomycin and maintained in a humidified atmosphere at 37•C and 5% CO 2 , either when initially cultured upon P24 multi-well culture plates for evaluating growth and death rates to adjust the model, or when cultured upon the microfluidic system prototype for final validation purposes. The medium was changed twice a week. Confluent cells were detached from the dishes by using Trypsin-EDTA (0.05% in HBSS, HyClone) when needed.
After counting with the support of a hemocytometer, both cell lines were seeded at 10.000 cells/cm 2 in a P24 multi-well plate to determine the cell proliferation rates for collecting information with the cell culture plates in similar conditions to those that would be applied upon the final validation in the organ-on-chip prototype, in which the same number of cells is placed in the different inlets. Imaging and counting was also performed with the support of a Leica DMIRB inverted microscope equipped with a digital Leica DC100 camera (Leica, Nussloch, Germany). Cell viability was monitored for 11 days in P24 multi-wells and in the organ-on-chip prototype, by calcein/propidium iodide staining to analyze the number of living and dead cells and adjust the cellular automata. Cell viability was determined using a calcein/propidium iodide dual-staining assay (Invitrogen, Molecular Probes). Briefly, the cell culture medium was removed and the cells were rinsed with phosphate-buffered saline. Next, 1 μM calcein and 2 μM propidium iodide were added in each well and incubated at 37 °C for 20 minutes. Fluorescence was evaluated using an inverted Leica DMIRB microscope equipped with a digital camera, Leica DC100 (Leica, Nussloch, Germany).

Results and Discussion
Simulation results and potential applications of the model. The procedure for defining the boundaries of cellular automata lattices directly from the CAD models is tested with different designs (see Section 2.3), upon which the collective dynamic behaviors of cells are simulated, to better expose the potential of this approach. Some results from the developed simulator performing within different devices are included in Fig. 4, which shows the results for: a) Cells seeded in opposite inlets and growing through a single-channel device; b) a detailed view of another simulation within the same single-channel device, in which two cells types -healthy (green) and invasive (blue)-are seeded and in which the addition of a drug (pink) and its diffusion, interacting just with the invasive cells, is also modeled; and c) different cells seeded in radially placed inlets and proliferating along the radial channels to meet in a central chamber. The healthy (green) and invasive (blue) cells, with different growth rates, can be appreciated and the incorporation of a drug (pink) and its diffusion is also modeled. The red points represent dead cells according to the defined death probabilities (P d = 0.3 in these trials). The cellular automata is programmed so as to ask the user to select the seed zones for placing the initial cells and cell types and the stepped dynamic process can be also modified to add, at a certain step (or time), a drug or reagent, which typically interacts with one of the different cell types. In our cases we are modelling the interaction between the drug (pink) and the invasive (blue) cells. The growth rates and the diffusion speed of cells and incorporated substances can be defined and adjusted according to experimental data (see Section 3.2).
A more realistic view, thanks to the use of Matlab's CAD import tools and its three-dimensional views, is presented in Fig. 5, in which the simulation results of collective cell behaviors within different organs-on-a-chips are presented. Figure 5a shows a 3D view of the interactions among different cell types along a single channel and Fig. 5b corresponds to a 3D view of the interactions among different cell types within a blood-brain barrier chip. To illustrate the whole modelling process, Fig. 6 presents a summary of the development of a simulation for a multi-chamber and multi-channel organ-on-a-chip device conceived for studying metastases. Figure 6a shows the computer-aided design (NX-8.5) of the system and Fig. 6b shows the CAD incorporated to Matlab, while Fig. 6c presents the related lattice with the allowed and prohibited zones. The final microfluidic system, used as preliminary validation prototype, is based on the simple design of Figs 4a and 5a (mainly a microchannel connecting two inlet wells) and is obtained by PDMS casting upon a laser stereolithography mold.

Preliminary model adjustment through cell culture experiments. Adherent cells (MC3T3
-osteoblast precursors-and N2A -neuronal phenotype-) are used to adjust the growth rates and obtain a preliminary validation, connected to actual experimental data, of the cellular automata working upon CAD-based lattices. The growth rates are calculated, as detailed below, by culturing cells in multi-well plates and monitoring cell proliferation, by taking detailed microscope images from the cultures at different times. A total of 10 fields (n = 3) for each cell type culture were photographed at days 1, 3, 5, 7 and 11. Final calcein propidium iodide staining let us analyze the number of living and dead cells and adjust the cellular automata by setting P d . Then, the calculated rates and probabilities are introduced in the simulator and applied to modelling the collective behavior of cells cultured within a simple microfluidic system, similar to that of Figs 4a and 5a, whose prototype is also seeded with cells to provide a experimental comparison between O-o-C simulation and O-o-C culture for a first validation of the simulator. These cell culture processes are detailed in Section 2.4 and their results further illustrated in Figs 7 and 8.
The information from the cell cultures is summarised in Fig. 7, which presents (in Fig. 7A-J) the evolution along 11 days of the MC3T3 (Fig. 7A-E) and N2A (Fig. 7F-J) cells in a representative position of the culture well, taken as example from the culture tests used to evaluate proliferation dynamics. All culture data are presented (in Fig. 7K-L) in the form of summary graph of the cell growth dynamics (Fig. 7K) and cell viability (Fig. 7L)  , it is possible to adjust the required iterations for the different cell types in accordance with the days under culture and, hence, adjust the simulator, by employing different step numbers to trigger the proliferation step (Fig. 2a) of each cell type. The proliferation steps for each cell type, which lead to simulated cell (2019) 9:14886 | https://doi.org/10.1038/s41598-019-51494-1 www.nature.com/scientificreports www.nature.com/scientificreports/ numbers according to experimental cell culture results and are finally used for simulating the O-o-C of Fig. 6, are included in Table 1. Subsequently, Fig. 8a presents the simulation results upon the microfluidic system using the adjusted growth rates and providing a simulated overview of the dynamic growth process along 11 days after seeding the cells through the microsystem inlets. In green, the N2A cells are presented, while the MC3T3 cells are shown in dark blue; besides, the red points represent dead cells. For comparative purposes, Fig. 8b shows the actual cell culture results upon a physical prototype of the multi-chamber organ-on-a-chip device, shown in Fig. 6d, again along a 11-day cell culture process. It can be appreciated that the proliferation of the MC3T3 cells in the microsystem is faster than that of the N2A cells especially at days 5 and 7 (probably due to a lower cell adhesion of the N2A cells to the material), as also predicted in the simulation (see Fig. 8).
Final discussion and future research proposals. Taking into account that cell culture materials and methods convey relevant cost, time and personal dedication investments, we consider that counting with simulators for predicting the collective behaviors of cells within lab-and organ-on-chip devices may prove interesting for the further expansion of the field.
We have shown that, once the cellular automata model is adjusted with the information obtained from conventional cell culture tests, the simulator helps to predict the way that cells will proliferate within the actual organ-on-a-chip system. Consequently, these simulations can prove useful for estimating the required materials, facilities and testing conditions (i.e. days under culture), so as to validate in vitro the real performance of Detailed view of another simulation within the same single-channel device, in which two cells types -healthy (green) and invasive (blue)-are seeded and in which the addition of a drug (pink) and its diffusion, interacting just with the invasive cells, is also modeled. Selected iteration when the drug is added and aspect after addition of the drug. (c) Different cells seeded in radially placed inlets and proliferating along the radial channels to meet in a central chamber. Healthy (green) and invasive (blue) cells with different growth rates can be appreciated and the incorporation of a drug (pink) and its diffusion is also modeled. Selected representative iterations are shown. In all cases, the red points represent dead cells according to the defined death probabilities. For each model different representative iterations are selected to illustrate the dynamic process. Each pixel corresponds to a single cell (typically 20 × 20 μm 2 ). Scale bars: 2 mm. (2019) 9:14886 | https://doi.org/10.1038/s41598-019-51494-1 www.nature.com/scientificreports www.nature.com/scientificreports/ innovative organ-on-chip devices and related microsystems. They can be used also for analyzing if a proposed design is adequate for studying cell-cell and cell-material interactions, in connection with disease modelling and depending on the available testing facilities. For example, the simulations from Fig. 8 helped us understand that, if we wanted to fully exploit the multi-chamber and multi-channel features of the proposed design with the cells under study, then we should extend the culture beyond the mentioned 10 days, which is challenging, or redesign with a finer or downscaled inlet, channel and chamber network. Design optimisation on the basis of the in silico studies is, therefore, also possible.
Among current limitations we can mention that the simulator works in a 2D environment with boundaries defined in accordance with the organ-on-chip designs and without considering cell migration, just proliferation of different species and death probabilities step-by-step. This proves enough when using adherent cells seeded at very low densities and for short culture times, but needs to be updated for including the possibility of migration and of a real 3D environment for more precise studies. We expect to solve the issue of the third dimension in a   quite straightforward way, by making use of common slicing software employed for additive manufacturing by digital light processing, which generate (layer-by-layer) black and white masks, similar as the ones used in this study, but starting from complex 3D CAD files without intermediate conversion to.jpeg format 30 . Approaches from other sectors may be also employed towards three-dimensional cellular automata 31,32 and synergies with other stochastic algorithms, such as the Monte Carlo method, may prove interesting for modelling additional phenomena within the organ-on-chip devices 33 .
Finally, we expect to enhance these simulators in the near future by using the results from simulations based on the finite-element method (FEM) as input for dynamically adjusting proliferation rates and probabilities within the lattice. In this way, fluid movement and shear rates within the microsystem and their effects on certain proliferation rates and differentiation phenomena may be modeled. The effects of thermal treatments (i.e. hyperthermia) or drug use and their impacts on cell survival may be also simulated by combining FEM and cellular automata, as we plan to study soon.

Conclusions
Counting with simulation resources for predicting the collective behavior of cells within lab-and organ-on-a-chip devices may support the design optimisation of these innovative and useful biomedical devices and in the experimental planning for their validation. Therefore, in this study we have focused on the establishment of a direct connection between the computer-aided designs of microfluidic systems, especially labs-and organs-on-chips, and the lattices for discrete cell modelling approaches aimed at the simulation of collective cell interactions, whose boundaries can be now defined directly from the CAD models. We have illustrated the proposal using a quite straightforward cellular automata model, applied it to simulating cells with different growth rates, within a selected set of microsystem designs, and validated our approach by tuning the growth rates of different cell types with the support of cell culture experiments and by checking the results with a real organ-on-a-chip system. We call the proposed procedure "game of life on a chip" and we envision that similar modelling approaches may support developers of these types of medical devices to optimise their engineering-design process and to support with experimental planning for their validation and a more straightforward regulatory compliance assessment by means of selected experiments and exhaustive simulations. Table 1. Proliferation steps for each cell type, which lead to simulated cell numbers in accordance with the experimental cell culture results and are used to fine-tune the simulator.