In silico characterization of cell–cell interactions using a cellular automata model of cell culture

Background Cell proliferation is a key characteristic of eukaryotic cells. During cell proliferation, cells interact with each other. In this study, we developed a cellular automata model to estimate cell–cell interactions using experimentally obtained images of cultured cells. Results We used four types of cells; HeLa cells, human osteosarcoma (HOS) cells, rat mesenchymal stem cells (MSCs), and rat smooth muscle A7r5 cells. These cells were cultured and stained daily. The obtained cell images were binarized and clipped into squares containing about 104 cells. These cells showed characteristic cell proliferation patterns. The growth curves of these cells were generated from the cell proliferation images and we determined the doubling time of these cells from the growth curves. We developed a simple cellular automata system with an easily accessible graphical user interface. This system has five variable parameters, namely, initial cell number, doubling time, motility, cell–cell adhesion, and cell–cell contact inhibition (of proliferation). Within these parameters, we obtained initial cell numbers and doubling times experimentally. We set the motility at a constant value because the effect of the parameter for our simulation was restricted. Therefore, we simulated cell proliferation behavior with cell–cell adhesion and cell–cell contact inhibition as variables. By comparing growth curves and proliferation cell images, we succeeded in determining the cell–cell interaction properties of each cell. Simulated HeLa and HOS cells exhibited low cell–cell adhesion and weak cell–cell contact inhibition. Simulated MSCs exhibited high cell–cell adhesion and positive cell–cell contact inhibition. Simulated A7r5 cells exhibited low cell–cell adhesion and strong cell–cell contact inhibition. These simulated results correlated with the experimental growth curves and proliferation images. Conclusions Our simulation approach is an easy method for evaluating the cell–cell interaction properties of cells.

endothelial cells show high intercellular adhesion and they ultimately form tight junctions [10,11]. Fibroblasts and smooth muscle cells adhere to each other and form cell aggregates [12,13]. In the process of epithelial-to-mesenchymal transformation, which is essential for animal and cancer development, intercellular adhesiveness decreases [14]. Increasing the density of endothelial or smooth muscle cells increased cell-cell contact and decreased cell spreading, leading to growth arrest [15]. Furthermore, the homotypic intercellular adhesive forces between ectoderm, mesoderm, and endoderm cells are different in each other [16]. Although these cell-cell interactions are important properties of cells, they are hard to evaluate in standard cell assay systems. Therefore, simple methods for evaluating cell-cell interaction properties are required for the development of cell assay systems for drug screening.
Mathematical modeling and computational simulation using differential equations and cellular automata assume an important role in experimental modeling to study a variety of biological events [17]. Modeling is performed by abstracting some characteristic factors from the biological events, and the simulation displays the reconstruction results when those factors are used. Thus, modeling and computational simulation can help us to understand how certain biological processes occur, as well as inform us of the critical factors involved in these biological processes. There are many cellular automata models that can describe the population dynamics of cultured cells [1,[18][19][20][21][22]. These models capture the features of the proliferation process and include important parameters (cell size, seeding density, spatial distribution of cells, migration, and oxygen density) necessary for prediction of cell population behavior. Furthermore, there are some cellular automata models that include cell-cell adhesion effects to simulate cultured cells in 2D and 3D [19,20,23]. Therefore, using these models, we can estimate the cell-cell interaction properties of cultured cells.
In this study, we developed a system for estimating the cell-cell interactions of cultured cells using a cellular automata model. To estimate the interaction parameters, we focused on images of cells in 2D cell culture because cell-cell interactions, especially cell-cell adhesion and cell-cell contact inhibition (of proliferation), affect the formation of cell aggregates. We obtained cultured cell images over time, and then assessed these images for corroboration with the cellular automata system. Using the cell growth rate and seeding density, which were determined from the images, we succeeded in estimating the cell-cell adhesion and cell-cell contact inhibition parameters. This method is useful for estimating these cell-cell interaction properties. This method could therefore be used as a new cell assay system for analysis of cell-cell interactions during drug screening.

Materials
Male Fisher 344 rats were purchased from Japan SLC (Shizuoka, Japan). The rat aorta smooth muscle cell line, A7r5, was obtained from DS Pharma Biomedical Co. Ltd (Osaka, Japan). The human cervical cancer cell line, HeLa, and the human osteosarcoma cell line, HOS, were obtained from the Health Science Research Resources Bank (Osaka, Japan). Cell culture medium was purchased from Sigma-Aldrich (St. Louis, MO). Fetal bovine serum (FBS) was purchased from JRH Biosciences (Lenexa, KS). Antibiotics were purchased from Life Technologies Japan Ltd. (Tokyo, Japan). Other reagents were purchased from Wako Pure Chemical Industries Ltd. (Osaka, Japan), Sigma-Aldrich, and Life Technologies Japan Ltd.

Preparation and culture of rat mesenchymal stem cells
Rat mesenchymal stem cells (MSCs) were isolated and primarily cultured as previously described [24]. Briefly, bone marrow cells were obtained from the femoral shafts of 7-week-old male Fisher 344 rats, which were anesthetized and euthanized by exposing of carbon dioxide. The cells were obtained from at least two rats and pooled in order to reduce the influence of individual differences. The culture medium was Eagle's minimal essential medium (with Earle's Salt and l-glutamine) containing 15% FBS and antibiotics (100 units/mL penicillin G, 100 µg/mL streptomycin sulfate, and 0.25 µg/ mL amphotericin B). The medium was replaced twice a week, and cells at passages 2-3 were used in this study. This study was carried out in strict accordance with the recommendations in the Guide for the Care and Use of Laboratory Animals of the University of Kitakyushu. The protocol was approved by the Committee on the Ethics of Animal Experiments of the University of Kitakyushu.

Cell culture
A7r5 cells, HeLa cells, and HOS cells were cultured in DMEM supplemented with 10% FBS and antibiotics (100 units/mL penicillin G, 100 µg/mL streptomycin sulfate). The medium was replaced twice a week.

Cell staining
Cells were seeded in a 35-mm culture dish at around 1 × 10 4 cells/cm 2 . The cells were fixed with 4% paraformaldehyde and stained with 0.4% trypan blue solution. The cells were imaged at 8.4× magnification using a stereomicroscope (SZX12; Olympus, Tokyo, Japan) equipped with a DP70 color charge-coupled device camera (Olympus).

Image extraction of cell distribution
The obtained cell images were analyzed using Image J software (NIH, Bethesda, MD). Each cell image was split into each RGB color channel. Then, the red channel image, which was the image with the highest contrast against the background, was subtracted from the background light shadow. The image was binarized with the adequate threshold intensity value, which was determined by referring to the highly magnified image. The size of the square image was dependent on each cell type (rat MSC, 3.5 × 3.5 mm; Hela and HOS cells, 3 × 3 mm; A7r5 cells, 4 × 4 mm). After clipping the square image from the binarized image, the resolution of the image decreased to 100 × 100 px. The clipping cell size was determined by the area of approximately 1 × 10 4 cells. We prepared more than five images for each condition.

Simulator specification
The 2D cell simulator was developed using Java (Oracle, Redwood Shores, CA). The cell distribution of the simulator was modeled by the cellular automata. The calculated space was set to 120 × 120 cell units and the displayed space on the screen was 100 × 100 cell units of the center of the calculated space ( Fig. 1). Each unit indicates whether it is occupied by cells (black point) or it is vacant (white point). The model includes cell movement and division regulated by cell-cell interactions (cell-cell adhesion inhibits cell movement and cell-cell contact inhibition inhibits cell division). One cycle of the calculation of cellular events (movement and division) indicates 10 min in a virtual environment. Because cells in the early stages of cell culture exhibit a proliferation lag, we simulated cell proliferation behavior after 24 h of culture. Therefore, the cell seeding number, s 0 , was used to represent the cell number after 1-day culture in each cell type.

Cell dynamics
Each cell unit contacts the neighboring 8 units (Fig. 2). According to a previous study using cellular automata dynamical simulation of cell culture [21], the influence probabilities for the center unit from the neighboring 8 units are 1/12 and 2/12 at a diagonal position (P inf× ) and a side position (P inf+ ) of a center unit, respectively (Fig. 2). Each influence probability was derived from the ratio of the central angle to 2π rad on a circle with a radius equal to the length of one side of the unit square. Figure 3 shows a flowchart of the cellular automata simulation. The method of cellular automata was applied in all cell occupied units.
The cell movement event is produced according to the probability, P mot , which depends on the cell motility parameter mot; P mot = 1/(6 mot). The mot means the time (h) for a single unit transfer of the cell on average. During the cell movement event, the cell can escape from the event according to the influence of the surrounding cells. The total influence of surrounding cells, P su , is determined based on the influence probabilities, a diagonal position of cell number n × , and a side position of cell number n + , P su = n × P inf× + n + P inf+ = (n × + 2n + )/12.
Using the cell-cell adhesion parameter a, the cell escapes from the movement event with the probability, P motesc , If the cell does not escape from the event, the cell moves around a vacant unit depending on the influence probabilities.
The cell division event is produced according to the probability, P div , which depends on the cell doubling time where α (=0.7147) is the offset of the cell division due to overlap with the immediately preceding cell division. During the cell division event, the cell can escape from the event according to the influence of the surrounding cells. The total influence is determined by the influence probabilities (P inf× , P inf+ ) and the cell number of the diagonal position and the side position of the cell, ci = 12 n × P inf× + 12 n + P inf+ = n × + 2 n + . According to the cell-cell contact inhibition parameters, the cell escapes from the division event with the probabilities P ci . The P ci can be fixed arbitrarily; but, in this study, we used four

Estimation of cell proliferation parameters
The simulated cell proliferation parameters could be estimated from the growth curve and cell proliferation images. Cell proliferation was serially simulated with various cell proliferation parameters. Then, the simulated cell number was evaluated and rated by comparing with the experimentally obtained data using least square analysis. Within the several higher conditions, the most matching parameters were finally determined by visually comparing the simulated cell images with the experimentally obtained cell images.

Analysis of cell proliferation under experimental culture conditions
First, we analyzed cell proliferation using four different types of cells: rat mesenchymal stem cells (MSCs), human cervical cancer HeLa cells, human osteosarcoma HOS cells, and rat aorta smooth muscle A7r5 cells. These cells were stained with trypan blue daily. The obtained cell images were binarized and clipped into squares of about 10 4 cells. Then, the resolution of the clipped image was decreased to 100 × 100 px. Figure 4a shows the binarized and clipped images of each cell type. The appearance of the proliferation process varied according to the cell type. The area covered by MSCs showed a mottled pattern at 4 days and cell-free area remained even at 7 days. HOS cells grew in large colonies initially, but eventually showed a confluent pattern. HeLa cells grew rapidly and quickly reached confluence. A7r5 cells showed sparse cell confluence after 7 days of culture. The growth curves, calculated using number of black pixels in the cell images, are shown in Fig. 4b. The cell-stained pixels grew exponentially in the early phase of cell culture. Therefore, the transition of the pixel number represents the cell proliferation curve. Next, we calculated the doubling time and the pixel number at 24 h for each cell type using the pixel proliferation data from the early phase of the cell culture. The calculated doubling times of MSCs, HOS, HeLa, and A7r5 cells were approximately 27.7, 19.7, 24.2, and 40.3 h, respectively. The pixel number at 24 h for each cell type was 1390, 740, 1450, and 1070, respectively. Using these experimental data, we analyzed the cell-cell interaction properties of each cell type.

Development of cell proliferation simulator
To estimate cell-cell interaction properties such as cell-cell adhesion and cell-cell contact inhibition, we constructed a simple cell proliferation simulator using a cellular automata model. Here, cell-cell adhesion and Fig. 2 Array of unit squares representing cells. One unit is in contact with the neighboring 8 units. l is the length of a unit, which is altered in each cell type. θ 1 (π/3) and θ 2 (π/6) are the angles which face to the side and the diagonal position of units, respectively. The influence probabilities from the neighboring units are 1/12 at a diagonal position (P inf× ) or 2/12 at a side position (P inf+ ), which is determined by the angle cell-cell contact inhibition show cell migration inhibition and cell division inhibition, respectively. The graphical user interface of our developed system is shown in Fig. 5. We arranged 120 × 120 cell units of virtual space as cell proliferation space (Fig. 1), and a single black point was used to represent a single cell (Fig. 2). We assumed that the centered 100 × 100 cell units of the virtual space reflected the binarized and clipped images of cultured cells (Fig. 4a). To construct the simulator, we set 5 parameters, which were initial cell seeding cell number s 0 , doubling time t d (h), motility mot (h), cell-cell adhesion a (%), and cell-cell contact inhibition P ci (detailed in "Methods").
First, we showed outputs from our simulator. The outputs of this simulator were the cell number and cell image for each culture time, based on setting the above five parameters. Figure 6 shows results when the initial cell number s 0 was variable. The growing cells gradually filled the cell culture space, eventually occupying the entire space (Fig. 6). Comparing low versus high s 0 values, the cells with low s 0 values grew and formed a large cell colony, whereas cells with high s 0 values completely filled the space (Fig. 6). Figure 7 shows results when the doubling time t d was variable. When t d was set to 18 h, the cells occupied almost the entire culture space on day 5 of culture. However, when t d was set to 48 h, the slope of the growth curve decreased and the cells could not fill the space by day 7 (Fig. 7). Figure 8 shows the results when the motility mot was variable. The growth curves demonstrate gradual growth when the mot value was high, i.e. low motility conditions (Fig. 8). Focusing on the growing cell images at day 4 of culture, the cells were well spread when the mot value was low, whereas they displayed a mottled pattern when the mot value was high (Fig. 8). Figure 9 shows the results when cellcell adhesion a was variable. The slope of the growth Fig. 3 Flow charts of the cell proliferation simulation. The left chart shows the main process and the right chart shows the cellular automata process. One cycle of the process is 10 min of virtual cell culture. There are two events, namely movement and division, in the cellular automata process. The movement event is started according to the P mot , which is determined by the cell motility parameter mot, and terminated according to the P motesc , which is determined by the effect of cell-cell adhesion (parameter a) with surrounding cells (the parameter is P su ). The cell division event is started according to the P div , which is determined by the cell doubling time t d , and terminated according to the P ci , which is determined by the effect of surrounding cells curve gradually decreased when the value of a was high (Fig. 9). The cell images obtained after 4 days of culture showed a highly mottled pattern when the value of a was high (Fig. 9). Figure 10 shows the results when cell-cell contact inhibition P ci was variable. When the P ci was set as null, the cells completely filled the space after 7 days of culture (Fig. 10). However, when the P ci was set to strongly inhibited conditions, the cells proliferated but only sparsely filled the culture space (Fig. 10).
The initial cell number s 0 and the doubling time t d are critical factors for cell proliferation, and when they are changed, the form and slope of the growth curve drastically changes (Figs. 6, 7). In this study, these two parameters were determined by the experimental results. However, the effects of motility mot, cell-cell adhesion a, and cell-cell contact inhibition P ci were relatively restricted. These parameters affected growing cell patterns (Figs. 8,9,10). In this study, we could not determine these parameters from the experimentally obtained results. Of these parameters, mot and a are parameters for assessing for the cell movement event (Fig. 3). Comparing the cell proliferation simulation with variable mot to that with variable a, the impact of a was more significant than that of mot (Figs. 8, 9). Furthermore, cell motility can be determined by some other methods [7,25]. Therefore, in this study, we simulated the cell proliferation behavior with the cell-cell adhesion a and cell-cell contact inhibition P ci , which are the parameters reflecting the cell-cell interaction activity, as variable.

Simulation of experimentally obtained cell proliferation
One of the advantages of simulation is obtaining hypothetical experimental data. We can serially simulate various cell proliferation behaviors with varied cell proliferation parameters. Therefore, we attempted to estimate the cell-cell interaction activity parameters, i.e. cell-cell adhesion a and cell-cell contact inhibition P ci , for MSCs, HOS, HeLa, and A7r5 cells. The experimentally obtained initial cell number and doubling time of each cell were used in the simulation. The motility was defined as a constant value. a varied from 10 to 40, in intervals of 10. P ci was defined as null, weak, positive, or strong inhibition. The simulated data were evaluated by comparison with experimentally obtained data. The simulated growth curves were ranked by least square analysis with the experimental growth curves. After choosing several highranking conditions, we visually compared the simulated cell images with experimentally obtained images, and then determined the simulation parameters of each cell. The growth curves of each cell type. The data shows the average pixel number ± standard deviation, which are calculated from more than five cell images for each condition Figure 11 shows the simulated cell proliferation behaviors of each cell type. The simulated growth curves fit with the experimental data for each cell type (Fig. 11). With regard to the cell growth images, simulated images roughly resembled the experimental images for each cell type. Simulated MSCs demonstrated a high cell-cell adhesion (a = 40) and positive cell-cell contact inhibition (Fig. 11a). These cells exhibited fine mottled patterns of cell growth at day 3, and they did not reach confluence at day 7 (Fig. 11a). The fine mottled pattern was similar to that observed in the experimentally obtained cell proliferation images. Simulated HOS cells showed low cell-cell adhesion (a = 20) and weak cell-cell contact inhibition (Fig. 11b). The simulated growth curve highly correlated with the experimental results. However, the experimental HOS cells exhibited a large colony cell pattern at day 3, but the simulation could not reproduce the cell patterns (Fig. 11b). The simulated HeLa cells demonstrated low cell-cell adhesion (a = 20) and weak cell-cell contact inhibition (Fig. 11c). The simulated growth curve and cell images correlated with the experimental results (Fig. 11c). The simulated A7r5 cells demonstrated low cell-cell adhesion (a = 20) and strong cell-cell contact inhibition (Fig. 11d). The simulated growth curve of A7r5 cells also correlated with experimentally obtained results. The simulated cell images showed sparse cell confluence with many vacant cell spaces at the later stages of cell culture, which was consistent with the experimentally obtained A7r5 cell pattern at day 7 (Fig. 11d).
In our system, the initial cell number and cell doubling time were experimentally obtained for each cell type. These two parameters most significantly affected the cell growth curve, and thus our simulated growth curves highly correlated with those that were experimentally obtained (Fig. 11). However, simulated images of growing cells were not always consistent with the experimentally obtained images. Cell-cell adhesion affected the fine mottled formation of cell patterns, and cell-cell contact inhibition affected cell confluency at a later stage of cell culture (Figs. 9, 10). Therefore, we successfully simulated MSCs, which formed fine mottled patterns; HeLa cells, which demonstrated global proliferation; and A7r5 cells, which reached highly sparse confluence (Fig. 11). However, we could not fully simulate HOS cell proliferation, which demonstrated a large colony growth pattern (Fig. 11). The formation of large colonies requires a small number of initial seeding cells (Fig. 6). However, the initial cell number of HOS cells was about 740 (Table 1), which is higher than the initial cell number forming the large cell colony (Fig. 6). This inconsistency is attributable to the cell attachment behavior to the culture substrate. The experimentally obtained cell images of HOS cells at day 1 showed many small cell aggregations that contain several cells (Fig. 4). The cell-seeding algorithm of our system was assigned a random number. Thus, the simulated initial cell location was more dispersed than the experimental HOS cell initial location. We then simulated the proliferation of HOS cells using the experimentally obtained cell image from day 1, and repeated the cell simulation with cell-cell adhesion a as variable and weak cell-cell contact inhibition. With low cell-cell adhesion (a = 10), the simulated growth curve correlated with the experimental growth curve of HOS cells (Fig. 12). Furthermore, the simulated cell image showed large cell colony forming proliferation (Fig. 12). Therefore, the cells that attached to the substrate with small cell aggregation were well simulated using the experimentally obtained initial cell image.
A summary of the simulated parameters for each cell type is shown in Table 1. From our simulation, we identified three modes of cell proliferation behavior based on the four cell types analyzed. The first mode includes The cancer cells demonstrated exponential proliferation and reached more than 95% confluence (Fig. 4). Generally, cancer cells can grow without any positive cell-cell adhesion or cell-cell contact inhibition, and they can easily attain an over-confluent cell pattern. Therefore, our simulated results reasonably describe these cells. The MSCs demonstrated a fine mottled proliferation pattern (Fig. 4), which requires high cell-cell adhesion values to be produced by the simulation (Fig. 9). Mesenchymal cells exhibit cell-cell contact inhibition and intercellular adhesive force is reported to be high [16]. Therefore, we concluded that the selected value of the cell-cell adhesion of MSCs was reasonable. A7r5 cells showed a very sparse confluence pattern (Fig. 4). These sparse patterns require strong cell-cell contact inhibition. On the other hand, simulated A7r5 cells did not demonstrate high cell-cell adhesion. This is likely because the growth curve of A7r5 cells was exponential and the proliferative cell pattern did not show a mottled pattern. Thus, we suppose that this is the interaction properties of A7r5 cells.
In this study, we simulated all cells with a constant motility (mot = 4.0). This is because the effect of motility was limited in the simulated growth curve (Fig. 8). It is difficult to determine appropriate parameter values when the number of variables increases. Cell simulated results, in particular, can be quantitatively compared to experimentally obtained cell growth curves, but only qualitatively compared to experimentally obtained cell proliferation patterns. Therefore, it is necessary to reduce the number of variables. In this study, we treated only two parameters, cell-cell adhesion and cell-cell contact inhibition, as variables, and we succeeded with cell proliferation simulations (Figs. 11, 12). Future studies on the quantitative comparison of simulated cell images and experimentally obtained images are needed. However, it is difficult to evaluate the similarity of two different images. Here, we attempted 2D Fourier transformation of cell images and compared the power spectrums of the resulting images, but we did not obtain quantitative values from the comparison. Therefore, another image comparison algorithm will be needed for this type of analysis.

Conclusions
We developed a cellular automata simulator for determining the cell-cell interaction properties of cultured cells. Cultured cells proliferated with their characteristic cell patterns and growth curves. We determined their doubling time from their growth curves, and then simulated their proliferation using their doubling times and initial cell numbers, with cell-cell interaction properties as variables. The proliferation could be simulated and we could determine their cell-cell interaction properties. HOS and HeLa cancer cells exhibited low cell-cell adhesion and weak cell-cell contact inhibition of proliferation. MSCs exhibited high cell-cell adhesion and positive cell-cell contact inhibition of proliferation. A7r5 cells exhibited low cell-cell adhesion and strong cell-cell contact inhibition of proliferation. These cell-cell interaction properties, which are characteristic features of each cell, could not be determined from typical cell culture experiments. Therefore, our developed simulation approach is an easy method for evaluating the cell-cell interaction properties of cells. Furthermore, this simulation can be used as a basic system for simulating many types of cellular events.   Fig. 12 Cell proliferation simulation of HOS cells using the experimentally obtained initial cell pattern. The upper graph shows the growth curves and the lower images are simulated (line and Sim) and experimentally obtained (circles and HOS) proliferating cell images. The doubling time t d was used as determined in Fig. 4. Other parameters were: mot, 4; a, 10; P ci , weak