In-vitro cell studies

The tumour-like hydrogel network model was successfully printed (Fig. 2a). To determine the proliferation of MDA-MB-231 cells embedded within the hydrogel, MTT assay was used to monitor cellular metabolic activity on days 0, 4, 7, 10, and 11. As illustrated in Fig. 2b, compared to day 0, MDA-MB-231 cells in the 3D cell/hydrogel construct demonstrated 1.86-, 2.7- and 2.78- and 2.8- fold proliferation on days 4, 7, 10, and 11, respectively. Indeed, cells showed rapid proliferation in the first 7 days and almost maintained almost a plateaued cell proliferation from day 7 to 11. Interestingly, the results showed that the doubling time of the MDA-MB-231 encapsulated in the 3D microenvironment was three times higher than that of cells grown in 2D cultures, which can be attributed to the reduced cell activities within the 3D matrix²². From day 7 to day 11, the number of proliferating cells remained approximately constant and left a question behind: whether cells died or entered to non-proliferating phase (quiescent) due to undesired conditions. To answer this question, live-dead assay, as well as Ki-67 immunostaining, were further employed during 11-days period to get a better understanding of cell behaviour growth encapsulated in the 3D scaffold.

(A): Fabricated cell/hydrogel structure using 3D bioprinting technique; (B): MDA-MB-231 proliferation embedded in hydrogel network post printing from day 0 to day 11. Error bars represent ± SD, n ≥ 3.

The viability of MDA-MB-231 over 11 days was visualized using live-dead staining, which was in line with MTT assay results. As shown in Fig. 2, most cells were viable after printing, and the viability rate of cells was 76 ± 2% on day 0, which clearly demonstrated minor damage of the bioprinting process on cell viability. The rate of viability increased over the first week and reached 98 ± 1% and 99 ± 1% on day 4 and day 7, respectively. Therefore, the structure was porous enough for oxygen and glucose to diffuse and distribute through the hydrogel scaffold, which could provide a proper living environment for cells. From day 7 to day 11, although some cells died, the majority of cells survived, and the rate of viability reached 96 ± 2% on day 11. By comparing the results of live-dead assay and MTT assay, it can be concluded that a significant portion of cells entered the resting phase after seven days since the maximum capacity of the scaffold had been achieved, and some of the cells began to die during the long-term stationary phase, or due to the lack of resources. Cell death in this experiment is almost negligible, which illustrates the promising potential of gelatin/alginate scaffold for bioprinting applications for a long time.

To visualize the proliferative capacity of MDA-MB-231, cells were fixed, and an anti-Ki-67 antibody was used to image proliferating cells using a confocal microscope (Figs. 3, 4). Ki-67 is a common-used marker that is present for all active phases of the cell cycle but absent in cells at the stationary phase²³. The results demonstrated that on days 0 and 4, almost 98 ± 1% and 95 ± 2% of the cells were positive for ki-67, respectively, while this number decreased to 86 ± 2% on day 7, followed by a dramatic drop to about 48.2 ± 2.4 on day 11. This result is in agreement with the data in the previous (Fig. 3), illustrating that within seven days, cells not only survive but also maintain their proliferating ability. From day 7 to day 11, although a high proportion of cells demonstrated to be alive, they were quiescent and were not able to proliferate anymore because of the lack of enough space for cells to proliferate. Additionally, on days 7 and 11, cells became more aggregated, particularly close to the pores, and cells at the center of the cell aggregates were shown to be non-proliferating due to being surrounded by other cells. Therefore, the proliferation process was aborted as there was not enough space for daughter cells to be placed and due to the lack of nutrients and oxygen at the center of aggregates.

Microscopic images demonstrating viability of MDA-MB-231 cells within 3D hydrogel constructs using the fluorescent live/dead assay kit from day 0 to day 11. Live and dead cells were stained using calcein-AM and PI, respectively (green color represents live cells; red color represent dead cells). Cells were imaged using the laser scanning confocal microscope. Scale bar, 50 μm (enlarged images, scale bar, 30 μm).

Ki-67 staining of encapsulated MDA-MB-231 cells within 3D bioprinted constructs. Cells were stained using anti-Ki-67 antibodies visualized with Alexa Fluor 546, and Hoechst 33,342 (red color represents cells positive to ki-67; blue color represents all cells). Cells were imaged using the laser scanning confocal microscope. Scale bar, 50 μm (enlarged image, scale bar, 30 μm).

The results of the in-vitro experiments were applied to parameterize the developed mathematical model.

In-silico cell studies

An essential property that makes the 3D bioprinting technique superior to the other 3D cell culture techniques is the precise control it offers over created structures. This property provides the opportunity to fabricate bio-mimetic structures with controlled structural and mechanical properties such as porosity, permeability and stiffness with high resolution^24,25,26. However, analyzing post-printing cell behaviour within the scaffold depends only on performing different in-vitro measurements. Experimentally measuring some of aspects of cellular behaviour in 3D bioprinted structures is challenging due to due to the lack of precise quantitative techniques, such as defining cell–cell and cell-microenvironment interactions. This problem has motivated our development of a mathematical framework to simulate post-printing cell behaviour within the scaffold. This framework must not only overcome these challenges but also accurately designs and predicts post-printing cellular functions without replicating experiments.

Individual-based modelling using a cellular automaton is one method to simulate the spatial and temporal mechanism of cell growth at the cellular level^27,28. This approach is a dynamic system that includes grids of cells, and each cell has sets of discrete states²⁹. CA modelling has been widely used in recent years to investigate different cancer cell mechanisms based on a variety of static automaton rules³⁰. However, to date, no study has applied CA modelling in 3D culturing cancer cells using 3D bioprinting. We selected this mathematical method for this study on simulating cell growth encapsulated in a 3D bioprinted structure due to its ability to capture spatial properties of the 3D printing structure and its flexibility to explore different hypotheses. Additionally, since the data obtained from our in-vitro experiments contained cells in discrete form, this discrete mathematical technique would have a more accurate simulation of this process. The framework developed in this study represents rules in cellular proliferation, viability, movement, and interactions with the environment, which includes hydrogel and neighbouring cells, as well as cluster formations within the hydrogel network. The in-silico results demonstrated in this article are based on mean values, and standard deviations from n = 100 simulation runs, where n is motivated by consistency analysis (Supplementary Material, Consistency Analysis).

Figure 5, showing the cell proliferation patterns for 11 days within the scaffold for both in-vitro and in-silico, illustrates them in agreement with each other. In cell proliferation, the initial cell density and scaffold capacity are two key parameters that we specified as \({C}_{\mathrm{initial}}\) and \(C\) variables, respectively. The corresponding values of these parameters were determined using calibration to produce the best fit to the in-vitro studies. Note that, the initial cell density in the developed model represents the effective initial population of cells that can interact with one another within a thin layer, not the total number of cells in the scaffold. Therefore, the simulated cell density was reduced by a scale factor compared to the cell density in the experimental settings. The results showed that cells reached the maximum cell density of cells in both simulation and experiments after seven days. The time it took for cells to achieve maximum density in the scaffold was strongly dependent on the number of initial cells and the capacity of the printed scaffold. The greater the number of initial cells and the less the scaffold capacity, the sooner cells reached maximum cell density and stopped proliferation. Therefore, by fine-tuning these parameters and running simulations in different scenarios, researchers can design the experiments to achieve desired results without repeating in-vitro assays.

Comparison between in-vitro and in-silico results of MDA-MB-231 cell proliferation within the 3D hydrogel network. Error bars represent ± SD, n ≥ 3.

Simulated data were also able to consistently replicate the viability and proliferation experimental results. As explained in the previous section, although cells showed around 99% viability within seven days of printing, the number of dead cells slightly increased from day 7 to day 11 when the significant portion of cells was in the resting phase. Hence, to precisely simulate the in-vitro condition, it was assumed that cells that remained in a prolonged stationary phase for more than the specified hours, defined by a stochastic number (\({C}_{d}\)), started to die with a specified probability (\({P}_{d}\)). This observation can be explained biologically by the cells’ inability to re-enter the cell cycle after entering the cell stationary phase.

Figure 5 shows the snapshots of in-silico MDA-MB-231 cells growing within the hydrogel network; as well as in-vitro microscopic images of cells in the 3D structure. In this figure, you can see cells’ distribution and progression of cell cluster formation on days 0, 4, 7 and 11 for both in vitro and in silico. In the in-silico images, yellow, red and black coloured cells are representative of proliferating, non-proliferating and dead cells, respectively.

On day 0, a few cells were distributed inside the hydrogel network. Over time, cells proliferated and created the first two-cell clusters and then bigger ones. Similar to in-vitro observations, the percentage of viability remained around 100% until day 11, when the viability slightly decreased and reached 93.74 ± 0.5%. Furthermore, simulated cell proliferation decreased over time and after seven days experienced a significant drop due to achieving the maximum capacity of the scaffold; and finally dwindled to 54.14 ± 0.25% on day 11. The animation of cell growth within the hydrogel scaffold in 11 days is also available in the Supplementary Material (Figure S3).

Another important factor apart from cell viability and proliferation is the ability of the cells to move in their surrounding matrix. This simulation can also be applied to analyze the cell movement as well as the structure and distribution of formed tumour clusters in the hydrogel network without experimental assessment. Tumour clusters might be created due to interactions between neighboring cells or between parent and daughter cells, depending on their position and the microenvironment³¹. Indeed, cells coordinate through cell–cell physical and signalling interactions and create clusters.

In Figure S1 (Supplementary Material), cells show a trend of crawling toward scaffold pores followed by forming clusters around those pores, as essential resources are in greater concentration there, particularly after seven days. This fact suggests that the hydrogel networks had limited resources transport capacity. Thus, we have defined particular ranges of attraction for both cell–cell signalling (\({L}_{C}\)) and cell-pore attraction (\({L}_{p}\)) to consider different cell migration directions. Speed of movement was another important parameter affecting cluster formation. Fallica et al. illustrated that the movement of cancer cells is inhibited in 3D microenvironments and shows extremely low speed due to the combination of material rigidity and the anchoring of cell receptors²². Therefore, based on the observations in this and previous studies, we have calibrated cell movement speed. In the movement processes, each individual moved at a specific time defined as \({m}_{C}\) in a direction determined based on cell attraction. During the process of calibration, by comparing in-vitro and in-silico results, it was concluded that cells had less tendency to move toward neighbouring cells within a smaller attraction range (\({L}_{C}\)) compared to the surface/pores of the scaffold (\({L}_{p}\)). Applying these rules in the model (Fig. 6), cells mimic the in-vitro cellular behaviour in terms of movement and cluster formation. The results of the proposed model are consistent with those of previous studies^22,32.

Middle panels visualize the MDA-MB-321 growth within the 3D hydrogel construct in silico; yellow represents proliferating cells; red represents non-proliferating cells; black represents dead cells. The right and left panels represent MDA-MB-231 cells encapsulated in 3D hydrogel constructs observed by a phase-contrast microscope on day 0, day 4, day 7, and day 11: scale bar, 50 μm.

In general, this model is developed to combine with in-vitro 3D-bioprinting evaluations, leading to a comprehensive analysis of the whole 3D fabricated structures. One of the main applications of this simulation is to predict the post-printing cellular behaviour in an unpracticed microenvironment which improves its capability to replicate desired biological settings. For example, this model provides the opportunity to evaluate the impact of different important parameters such as various initial cell densities on cellular behaviour in a long-term period. This can be of benefit to researchers to generate a more suitable microenvironment for cell growth without the need to repeat experiments. For instance, they can design the scaffold in terms of size or structural shape with the purpose of modifying the scaffold capacity to improve cell proliferation and decrease cell death.

In-silico model validation

To further validate the in-silico model, we performed the bioprinting procedures with different experimental variables in two different situations: case 1: varying initial cell densities; case 2: varying bioink formulation. In case 1, we did bioprinting with a bioink with 4% (\(w/v\)) gelatin, 4% (\(w/v\)) alginate, and 1.5 × 1 \({0}^{6}\) MDA-MB-231 cells \({\mathrm{mL}}^{-1}\). In case 2: we performed bioprinting experiments with a bioink with final concentrations of 4% (\(w/v\)) gelatin, 5% (\(w/v\)) alginate, and 2 × 1 \({0}^{6}\) MDA-MB-231 cells \({\mathrm{mL}}^{-1}\). Using the calibrated in-silico model, we would like to predict the proliferation pattern of cells in these two new conditions.

Figure 7, showing the cell proliferation patterns for 10 days for case 1 for both in-vitro and in-silico, illustrates them in agreement with each other. In the in-silico model, all parameters except \({C}_{\mathrm{initial}}\) (the initial cell density) have the same values as in the calibrated model. The simulated cell density is reduced by the same scale factor compared to the cell density in the experimental settings and set to \({C}_{\mathrm{initial}}=2000\). Both simulations and experiments demonstrated that cells did not achieve the maximal cell density after 7 days and kept growing. Hence, when the number of beginning cells reduced, the later cells attained their scaffold capacity and consequently ceased proliferating.

Prediction of cell proliferation using in-silico model for case 1: 1.5 × 1 \({0}^{6}\) MDA-MB-231 cells \({mL}^{-1}\) in 4% gelatin/4% alginate bioink; on day 0, 3, 7, and 10. In vitro study was performed using the MTT assay. Error bars represent ± SD, n ≥ 3.

Figure 8, comparing in-vitro and in-silico cell proliferation patterns for case 2, also shows agreement. In this case, experimentally, we altered the formulation of bioink. Increasing the alginate concentration can increase the rigidity of a hydrogel-based construction, as demonstrated earlier³³. It has also been found that the stiffness of the microenvironment would also affect cell movement and spheroid formation within the scaffold³⁴. Although parameters directly related to stiffness have not yet been integrated into our model, we may regulate cellular behavior and investigate the impacts of bioink formulation and stiffness on proliferation and migration by varying some rules of cellular movement. In the calibrated model for bioink containing 4% (w/v) gelatin and 4% (w/v) alginate, cells move every 15 h, denoted by \({m}_{c}\). Thus, with 4% (w/v) gelatin, 5% (w/v) alginate-based bioink, we reduce the movement speed of cells encapsulated in a stiffer microenvironment and change \({m}_{c}\) to 20 h while keeping other parameters unchanged. Comparing the results of an in-silico model to in-vitro data, we conclude that for 4% (w/v) gelatin and 5% (w/v) alginate-based bioink, \({m}_{c}\)=20 h closely match the in vitro proliferation trend within 11 days.

Prediction of cell proliferation using in-silico model for case 2: 2 × 1 \({0}^{6}\) cells \({mL}^{-1}\) in 4% gelatin/5% alginate bioink. In vitro study was performed using the MTT assay. Error bars represent ± SD, n ≥ 3.

The in vitro observations revealed that on day 11, cell proliferation decreased slightly, which can be explained by the stiffer microenvironment. Indeed, the rigidity might reduce cell movement and proliferation; and hinder the transport of nutrients, resulting in cell death over time. For more significant changes in bioink formulation, it is necessary to incorporate microenvironment stiffness or bioink-related parameters into the model to anticipate cell behaviour accurately. However, with 4% (\(w/v\)) gelatin, 5% (\(w/v\)) alginate, and minor modifications to the bioink formulation, our developed model can be used successfully.

Taken together, we could validate our model by creating variations in the in vitro data and successfully simulating the varied situation. Hence, we can confidently claim that this model can help researchers plan experiments more accurately by predicting the outcome. In fact, researchers executing simulations under various situations and fine-tuning related parameters can design experiments to reach the desired results without repeating in-vitro procedures.

Prospect

The developed mathematical framework in this study can be extensively applied in different bioprinting-related studies for various applications such as tissue engineering, oncology, and the pharmaceutical industry by extending its rules and improving its ability to provide an accurate prediction of biological systems. Our model can be expanded by incorporating bioink-related parameters such as stiffness and structural integrity, which regulate cellular behaviour, including proliferation and migration and oxygen/nutrient diffusion to the 3D network^{35,36,37,38,39}.

Moreover, the proposed model can be integrated with the machine learning algorithms and provide researchers with this opportunity to predict the temporal or structural effect of the hydrogel network on any desired objectives in the biological system. Furthermore, we can use CA simulation to pre-train the ML algorithm and then a transfer learning approach can be applied to train the experimental data.

Another prospect of this model is its application in a heterogeneous environment with multiple cell lines for studying cell–cell interaction and cell-ECM interactions. Besides, by adjusting the rules, this model can be integrated with pharmacokinetic modelling techniques to simulate drug treatment responses in 3D cell culture using 3D bioprinting to help study tumour development and metastasis, drug screening, and other aspects of cancer research.

In the end, we believe that this work or its combination with other modelling techniques can significantly influence the development of 3D bioprinting in the future and avoid conducting costly and time-consuming experiments to a great extent.