Abstract
Articular cartilage is characterized by low cell density of only one cell type, chondrocytes, and has limited self-healing properties. When articular cartilage is affected by traumatic injuries, a therapeutic strategy such as autologous chondrocyte implantation is usually proposed for its treatment. This approach requires in vitro chondrocyte expansion to yield high cell number for cell transplantation. To improve the efficiency of this procedure, it is necessary to assess cell dynamics such as migration, proliferation and cell death during culture. Computational models such as cellular automata can be used to simulate cell dynamics in order to enhance the result of cell culture procedures. This methodology has been implemented for several cell types; however, an experimental validation is required for each one. For this reason, in this research a cellular automata model, based on random-walk theory, was devised in order to predict articular chondrocyte behavior in monolayer culture during cell expansion. Results demonstrated that the cellular automata model corresponded to cell dynamics and computed-accurate quantitative results. Moreover, it was possible to observe that cell dynamics depend on weighted probabilities derived from experimental data and cell behavior varies according to the cell culture period. Thus, depending on whether cells were just seeded or proliferated exponentially, culture time probabilities differed in percentages in the CA model. Furthermore, in the experimental assessment a decreased chondrocyte proliferation was observed along with increased passage number. This approach is expected to having other uses as in enhancing articular cartilage therapies based on tissue engineering and regenerative medicine.
Introduction
Articular cartilage is a specialized connective tissue located at the ends of bones in the articular surfaces [1]. The functional properties of articular cartilage are related to reducing friction coefficient and support mechanical loads between opposing joint surfaces [2]. Given its damping function, articular cartilage is continuously exposed to mechanical loads that often result in an important number of individuals affected worldwide with cartilage related diseases [3]. Limited articular cartilage self-healing properties entail novel therapeutic options for its treatment, particularly when it is affected by traumatic injuries [4–7]. Among these are non-invasive interventions, such as physical therapy to relieve pain locally [8] or invasive treatments such as microfractures, mosaicplasty, osteotomy and Autologous Chondrocyte Implantation (ACI) to repair or reconstruct the injured cartilage [3, 9]. The ACI is a technique intended to expand the patient’s chondrocytes through in vitro cultures. Once a sufficient number of chondrocytes have been obtained, cells are re-implanted into the injured articular cartilage for restorative benefits [5, 11]. The employment of the ACI technique entails the consideration of several variables in the cell cultures such as the cell culture time, the cell population to be reached and the molecular synthesis [12]. Therefore, this treatment could be enhanced if computational methods allow to predict chondrocyte cell dynamics: migration, proliferation and cell death in order to improve in vitro procedures and increase the success of the therapy on articular cartilage recovery. The combination of agent-based models with experimental procedures, such as the use of a cellular automata (CA), can help to understand the behavior of experimental procedures to improve the cost-efficient performance of the experimental procedures carried out in cell culture [13].
A CA is an agent-based model that considers particles occupying a space in a regular lattice. Particles move randomly and change of state depending on current position and neighbors [14]. The CA can simulate cell dynamics such as migration, proliferation and cell death, provided that the CA is combined with stochastic tools that homogeneously model and distribute the particles in a given space [13]. In this context, the random-walk theory is a stochastic process formed by a successive sum of variables, independent and identical, randomly distributed (Gaussian distribution) [15]. Such theoretical approach has facilitated the simulation of the random cell migration, the number of proliferating cells during cell dispersion and the differentiation process [13, 16]. Moreover, the random-walk theory presents advantages over other models, such as diffusion. This model, used to simulate migration and proliferation, tend to create a smooth variation in cell density, which is an event that does not occur physiologically [16]. Additionally, the cell culture dynamics such as cell culture density, population-average, speed of locomotion and spatial cell distribution, all of them essential to determine temporal evolution of cell proliferation rates, can be predicted using a CA [17].
Regarding the applications, a CA helped to demonstrate that cell migration is a random fluctuation that does not depend on the surface area, but it is highly dependent on the initial parameters [18]. Longo and co-workers implemented a CA to analyze the spatial and temporal movements of blastocoel roof cells during epiboly. Their results showed that the model was able to predict cell rearrangements during epiboly and lateral dispersion of a patch of cells implanted in the blastocoel roof [19]. On the other hand, a CA in combination with random-walk theory was implemented to simulate the dynamics of muscle satellite cells: migration, proliferation and differentiation. Results from this work demonstrated that it was possible to predict the cell behavior through the computational model and the feedback with the in vitro experiments [13]. In this context, in order to validate CA results it is necessary to perform in vitro experiments for each cell type [13, 21]. Even though computer models have assessed cell dynamics in different cellular types [13, 21], to our knowledge, a CA model is yet to be done to simulate chondrocyte behavior. For instance, a CA was implemented to model the random cell motion and the production of substrate adhesion molecules of avian limb bud precartilage cells. Results showed that the CA adequately simulates the cell proliferation, the cell focal condensation over the well plate surface and substrate adhesion molecules production [22]. Additionally, the model was enhanced by adding multipixel cells with variable geometry and separation of spatial scales in the form of modeling cells and extracellular matrix molecules on grids of differing mesh size [23]. In these works, the implemented CA was used to simulate the chondrogenesis process, where mesenchymal cells proliferate and differentiate into chondrocytes forming a cartilaginous mold that subsequently undergoes in an ossification process (endochondral ossification) [24]. However, no studies have been carried out to assess the behavior of articular cartilage cells cultured in vitro. Thus, the aim of this research was to shed light on articular chondrocyte monolayer culture behavior. Accordingly, in this research the cell dynamics of chondrocytes was simulated using a CA, based on random-walk theory, with the aim to predict migration, proliferation and cell death. We hypothesize that the three processes can be accurately modeled as being completely stochastic.
Additionally, in this study two cell culture passages were assessed in order to document the replicative senescence experienced by a monolayer cell culture. Cell culture passages are referring to reseeding a subset of the cells obtained from previous proliferation periods in a new cell culture well plate. Furthermore, the surface tension effect was monitored in order to ensure that cell distribution was completely homogeneous over the entire surface of the well plate, since cells cultured in monolayer tend to locate in the well plate edges [13]. In this context, the results obtained from this study are relevant to enhance culture conditions of articular chondrocytes for future therapeutic strategies in tissue engineering and regenerative medicine. Moreover, the experimental and computational data allow to establish the best cell culture time and cell culture passage in order to standardize the in vitro procedures [6].
Material and methods
Experimental data
Cell culture
Human articular chondrocytes of the knee (LONZA, Walkersville, MD USA) were cultured and expanded at passages one (P1) and two (P2). Cells were then seeded onto 96 well plates at 3.000 cells per well (according to instructions by LONZA, Walkersville, MD USA). Culture media was changed every other day, and cells were maintained in humidified atmosphere at 37°C with 5% CO2.
Surface tension assay
Different culture media volumes per well were added in order to assess the surface tension effect on cell distribution within the well plate (Fig. 1A–E). Photographic records were daily captured within the week of culture. The meniscus formed at the superior border of the cell culture medium is the result of the surface tension effect. In this study, instead of measuring the surface tension, we rely on the theory of this physical phenomenon to explain the meniscus effect over the cell distribution.
Representative scheme to evaluate the surface tension effect over cell distribution. Different volumes of cell culture medium were implemented per each well plate. A) 50μl. B)100μl. C) 150μl. D) 200μl. E) 250μl.
Cell migration was evaluated by photographic record in the third day, when cells are in a 50% of confluence and it is possible to observe random movement. The chondrocytes were assessed under normal cell culture conditions (200μl of cell culture medium). Photographs were taken every two hours for a total period of 12 hours and the snap shot location was ensured by a reference point in the well. Pictures taken at different times can only illustrate cell movement inside the well plate. There is no quantification of how many cells and which distance they moved.
Proliferation and cell death assay
Cell proliferation was assayed on a daily basis for one week. Cells were counted with a hemocytometer; viability and cell death were determined by Trypan blue stain. Viable cells selectively inhibit Trypan blue to cross the cell membrane. Hence, dead cells become blue and can be counted under a light microscope. Proliferation between P1 and P2 was compared by slope values during the Lag phase: from day zero to day four of culture and log phase from day four to one week of culture.
Statistical analysis
All analyses were performed as mean ± SD (n = 3). Experimental and computational data obtained were analyzed by non-parametric statistics for each time-point (day zero, four days, one week). Wilcoxon rank test was performed for cell migration and proliferation. Kruskal-Wallis test for cell death (GraphPad v.6.0, Software Inc., California USA). In addition, in order to compare the proliferation differences between P1 and P2 chondrocytes, an F test was performed between slope values. Statistical differences were considered significant at p < 0.05.
Computational implementation
Cellular automata algorithm
To solve the CA the following algorithm was implemented (Fig. 2). First, the geometrical model (surface area and seeded cells) was represented and plotted in a bi-dimensional space (x-y plane). Secondly, an initial population of 3.000 cells was randomly distributed in 384.400 elements. Each element was only occupied by one cell. Once cells were arbitrarily distributed, a scheduled routine localized possible neighbors for each cell. Last, the model considered a computational framework to weigh probabilities for cell migration (p m ), proliferation (p p ) or death (p d ). The cells were randomly chosen.
Algorithm used to solve the CA. The algorithm consist in define the geometric area (well plate) and the initial cell population. Then, the cells are randomly placed in order to initiate a selective process for the cells can migrate, proliferate or die. If any of this events is executed a new cell population is generated, if not the cell remain in the same place. This process is repeated iteratively until complete one week of simulation with the aim of obtains a final cell population.
Due to the computational model included random cell movement, proliferation and cell death using a classical algorithm, the cellular dynamics was defined by weighted probabilities derived from experimental data using as the input parameter cell proliferation; therefore, the cell dynamics varies according to the cell culture period. According to this, the cell culture periods assessed in this study are the Lag phase (from day zero to day fourth) in which cells becoming accustomed to the ambient conditions, and the Log phase (from day fourth to one week) in which cells experienced an exponential growth. Therefore, depending on whether cells were just seeded or proliferated exponentially, culture time probabilities differed in percentages in the CA model. The outputs of the computational model are the weighted probabilities, expressed as a percentage, of the number of cells that are migrating, proliferating and dying. The computational analysis was carried out for both P1 and P2. The CA algorithm was performed using FORTRAN 90 (Formula Translating System, New York, USA) and results were visualized in TECPLOT 360 (Tecplot Inc. Bellevue), respectively.
The surface area of 96 well plates was represented as a circle of 0.3 cm2 (Fig. 3A). This area was gridded with a lattice of 10×10μm since it was assumed that a seeded chondrocyte should occupy this specific space according to cell dimensions described in the literature (Fig. 3B) [25]. Cell morphology was assumed as a square inside the 10×10μm grid for its representation in the mesh (Fig. 3C).
Geometrical model. A) A well plate from a 96 cell culture well plate with a superficial area of 0.3 cm2. B) A mesh created with Cubit 10.0 representing the elements where the cells are cultured. C) Mesh zoom showing the area of each cultured cell, called element.
According to the computational model, cells could migrate, proliferate or die. We call these cell dynamics events as cellular state. These states differed during Lag or Log phase. We attributed weighted probabilities as distinct percentages for each cellular state in order to fulfill the condition ∑p = 1, where p depends on migration probability (p m ), proliferation probability (p p ), and death probability (p d ), resulting in ∑p = p m + p p + p d = 1. According to Lag or Log phase, different probabilities were weighted for each cellular state.
Neighborhood characterization
To implement the CA model it is necessary to define a neighborhood around the element. Neighborhoods delimit the adjoining set of cells and relative position to each of them as a spatial arrangement of elements, where every single participant has an influence on one or more elements. Starting from an initial state of the cell population, cells change their states synchronously at each instant. This change is established via Von Neumann neighborhood system. The probabilistic aspects of the simulation must be considered within its dynamics.
Cell migration (p m , m d , a s )
Cell migration was modeled based on random-walk theory [13, 16]. Two parameters were additionally taken into account: cell migration distance (m d ) and available space (a s ). In the CA model, m d = 10 μm according to the previously described element dimensions [25]. If cell migration occurs, the CA model will consider the possible elements available to be occupied (Fig. 4A). In this case, a s = 1, 2, …, 5, since the cell may migrate onto four different elements or stay in the same place. According to previously established weighted probabilities, p m was determined by means of the condition . If all surrounding positions were available, (p m ) correspond to 1/5. On the other hand, if only one free element was available, the probability to fill the element would be one.
Cell proliferation, migration and cell death according to random-walk theory. A) Possible migration sites depending on probability and the neighborhood system of Von Neumann. A cell can migrate in all four directions (above, right, below or left) or stay in the same place. Due to optimization of the CA model cells cannot migrate in diagonal. B) Possible daughter cell states after proliferation depending on probability. A cell can proliferate in all four directions (above, right, below or left) or stay in the same place. Due to optimization of the CA model cells cannot proliferate in diagonal. C) Representative sample of cell death depending on probability. A cell can disappear leaving the space available (area in circle).
Cell proliferation was also modeled using the random-walk theory [13, 16]. This theory proposes a stochastic process, in which a cell can occupy a particular element according to the cell proliferation probability p p (Fig. 4B). In order to fulfill the requirement vacant positions were given by n = 1, 2, …, 5. Following the same approach as migration, if all neighboring positions were available the probability that a cell could proliferate and occupy an adjacent element was1/5. Furthermore, if only one element is available, the probability will be equal to one. Moreover, if all elements are occupied, it will represent a non-proliferating cell. As described for migration, weighted probabilities for p p depended on cellular phase: Lag or Log.
Cell Death (p d )
Finally, cell death was modeled based on the remaining weighted percentage depending on Lag or Log phases respectively. If a cell did not migrate or proliferate, it is determined that the cell would die p d = 1 or remain in the same element p d = 0 (Fig. 4C).
Results
Surface tension
Results of surface tension effect illustrate that culture medium volume had an effect on cell location (Fig. 5). It can be observed that cells at day one of culture have a homogeneous distribution over the well plate surface; however, after one week of culture the cells have a different behavior. For instance, at the lowest volume (50μL) cells tend to localize at the edge of each well in the plate (Fig. 5A); while at greater volumes (200μL), cell distributed homogeneously throughout the surface area of the well plate (Fig. 5B). In this study there was not a quantification of cells that are distributed at the edge of each well in the plate.
Surface tension effects on cell distribution. A) Cell distribution with 50μL of culture media. B) Cell distribution with 200μL of culture media. Scale bar = 500μm. Photograph at 4X magnification.
Experimental results showed random chondrocyte movements over the entire surface of the well plate (Fig. 6 A–C). Moreover, migration by weighted probabilities was derived from remaining computational data (Table 1). Cells cultured at P1 during the Lag phase migrated at a similar percentage compared to cells cultured at P2. On the contrary, cells cultured at P1 during the Log phase decreased migration percentages compared to the Lag phase. Cells cultured at P2 retained a similar migration percentage during the week of culture. The percentages shown in Table 1 represent the number of cells that are dedicated to migration. This means that the remainder percentage, either of the Lag Phase or Log phase, is the value needed to complete the 100% of cultured cells. Thus, the remainder percentage corresponds to the cells that are proliferating or dying. According to the percentages shown in Table 1, cells cultured at P1 tend to decrease their migration activity from 78% to 63%. This behavior indicates that cells are increasing their proliferative activity; therefore, there are fewer spaces for the cell to migrate. On the contrary, cells cultured at P2 tend to maintain their migration activity during Lag and Log phases (from 79% to 78%). It means that the amount of available spaces tend to be the same during the cell culture period.
Cell migration for cell culture passages. Representative sample of cell location. A) Cell location at time zero (8:00 am). B) Cell location after two hours (10:00 am). C) Cell location after four hours (12:00 m). Colored circles allow to evidence cell movement in each time period. Scale bar = 100μm. Photograph at 10X magnification.
Cell migration state. Weighted probability for P1 and P2
No significant differences were observed for all time points comparing the experimental data with the computational results (Fig. 7). On one hand, one difference between the experimental and computational results were found in cells cultured at P2, especially at day 2 of culture, which shows an increase in cell population; however, at this point, no significant differences were obtained. On the other hand, it is possible to observe the difference in the final cell population of the two cell culture passages. For instance, the final number of cells, for both the experimental and computational data at P1 was about 40.000 cells (Fig. 7A); whereas for cell cultures at P2 the final cell population were approximately 20.000 cells (Fig. 7B). This observation can be corroborated with the Table 2, where it is possible to observe the significant differences between the cell proliferation per day at P1 and P2. In addition, the cell distribution for both the cell culture and the computational model revealed a homogenous allocation (Fig. 8).
Cell proliferation curves for cell culture passages and CA. Representative sample for cell density. A) Cell proliferation for P1 vs CA. B) Cell proliferation for P2 vs CA.
Proliferation behavior per hour between P1 and P2
Proliferation behavior between chondrocytes culture atP1 and P2 for Lag phase was determined by calculating slopes values between day 0 and day 4 of culture *p = 0.004. For Log phase slope value between day 4 and one week of culture was highly significant ***p = 0.0001.
Cell proliferation comparisons between cell culture and computational model. A) Representative sample of cell distribution at P1. B) Representative sample of cell distribution at P2. Scale bar = 100μm. Photograph at 10X magnification.
Regarding the computational results, the proliferation rate was determined for each time point of cell culture. In this context, from the day zero up to 72 h of culture, similar percentages were observed for cell proliferation. In this time point the CA model predicted a 23% of proliferating cells at P1; while an 18% of the cells proliferated at P2. However, at the end of the culture (one week) cellular state differed with 36% of the cells proliferating at P1, and 21% proliferating at P2. This observation validates the use of this computational approach to accurately predict the cell proliferation of a given population.
No significant differences were observed among time points (Fig. 9). Experimental cell death was null on day zero when the cells were seeded. By the fourth day of culture cell death increased to 8.7%, and after one week of culture it was null again. Computational percentage defined a decreasing trend for cell death during culture (Table 3). In the Lag phase both cell cultures at P1 and P2 decrease from 5% to 2%; while in the Log phase both passages have a cell death percentage equal to 1%.
In vitro cell death determination. Cell death was calculated for cells at P2.
Cell death state. Weighted probability for P1 and P2
The aim of this research was to implement a CA model based on random-walk theory to predict the cell dynamics of chondrocytes: migration, proliferation and cell death. An accurate prediction of cell dynamics was possible thanks to the constant feedback between in vitro procedures, which established the input parameters – that indicate the amount of cells that are proliferating–, and the CA modeling. Accordingly, the cell proliferation allowed to predict the weighted probabilities of the number cells that are migrating, proliferating and dying at different culture times. In this context, the implemented computational approach can be extrapolated to enhance therapeutic techniques such as the ACI, due to the fact that CA can predict the final cell population in a given culture period. Thus, it might be possible to use ACI to approximate the in vitro experimental culture procedures, reducing costs and standardizing procedures [12, 26].
Regarding the cell dynamics, the photographic record showed cell movement into the well plate; however, as a study to quantify the amount of moving cells was not performed, the random movement of the cells was corroborated with [25]. In the CA proposed, a sensitivity analysis of cell migration distance (m d ) and available space (a s ) was done in order to elucidate the contact inhibition process to achieve better control of the migration and proliferation process. In this context, if the cells migrate and the neighboring elements are available, there will be a high number of proliferating cells; whereas, if the cells migrate around their original position, the quantity of proliferating cell will decrease. Therefore, m d and a s are fundamental parameters to determine the number of cells and homogeneous distribution inside the well plate [13]. In fact, it is not possible to compare the in vitro migration with the in vivo behavior, since chondrocytes in articular cartilage are surrounded by a specialized extracellular matrix which generates pressure on the cell and restrict its movement [25]. On the other hand, the geometrical form of the cell due to the movement into the well plate was not studied in this work, as the CA could represent only the state in a point. However, experimental results evidence the formation of typical lamellipodia: external protuberances on the mobile edge of the cell, when cells are migrating over flat surfaces (Fig. 6). This observation is consistent with results obtained by other authors [25]. The appearance of these protuberances (lamellipodia) is caused by the morphological change of the chondrocytes in monolayer cultures. Thus, chondrocytes change its round phenotype and tend to be similar to a fibroblast [27].
In this work was possible to compare cells at two different passages, observing a different cell dynamics. As it was mentioned above, chondrocytes cultured at P1 had a significantly higher proliferation rate compared with cells cultured at P2. This observation is consistent with those reported in the literature indicating that in vitro cultures experience replicative senescence even from P1 [28]. This has an important therapeutic implication, since cell passage or population doubling should be considered when implementing cell therapy [29].
In addition to cell dynamics, in this study was evaluated the surface tension effect on cell distribution. This variable is avoidable, unless cells are seeded at very low culture medium volumes. Additionally, cells cultured in a low culture medium volume, seem to be influenced in the mitotic response of chondrocytes. This observation can be corroborated with Fig. 5A, where it is possible to observe a less populated area in the middle of the well plate after one week of culture. Since surface tension effect was not modeled computationally in this study, to doing this it would be necessary to simulate the case when cells are located at the edge of the well plate. This phenomenon was not taken into account since the aim of this study was focused on normal culture conditions.
Even though more sophisticated models have been devised, such as Cellular Potts Model (CPM) that takes into account a more realistic framework including cellular morphology [30], this level of sophistication was not necessary to evaluate chondrocyte cell dynamics. As it was mentioned above, ACI requires a cellular expansion in a monolayer culture. For this reason, the CA that was implemented in this research works under a 2D environment and efficiently predicts chondrocyte dynamics. However, the methodology carried out in this work to develop the computational model can be extrapolated to perform models and imitate 3D tissue formation. For instance, some researchers have used the random-walk theory to simulate the random movement of cells into 3D models, predicting tissue growth [31, 32]. According to this, each computational model, either 2D or 3D, can use the same methodology and obtain different outcomes. Moreover, in order to validate the computational model, an exclusive experimental design must be implemented in each case. In this context, the present computational model represents an interesting bi-dimensional tool for predicting individual cell behavior in a normal monolayer culture.
Conclusion
This research evidences that a classical model, hand in hand with an experimental feedback, is sufficient to help predict the chondrocyte cell dynamics. In this context, the computational modeling using a CA is a promising tool that can be used to understand the behavior of chondrocyte culture procedures applied to articular cartilage repair. Moreover, the analysis of in vitro cell behavior is crucial for basic research as well as for quality control in cellular therapy [29]. Therefore, the procedures carried out in this research are suitable to be extrapolated to study cell behavior in treatments such as the ACI, representing new developments in tissue and regenerative engineering research.
Footnotes
Acknowledgments
The authors gratefully thank the research support group from the Republic Bank of Colombia – Foundation for the Promotion of Research and Technology (FPIT), professor Carmen Alicia Cardozo de Martinez and professor Afife Mrad de Osorio at the Institute of Biotechnology, Universidad Nacional de Colombia, for providing the lab space at the Biomimetics laboratory and the reactants to perform the experimental approach of this study.