1. Introduction
Non-small cell lung cancer (NSCLC) is one of the deadliest forms of cancer worldwide. Lung cancer organoids are three-dimensional cell aggregates grown from tumor cells in vitro. Unlike traditional monolayered cell cultures, patient-derived tumor organoids (PDTOs) preserve cellular architecture, mutations, and growth of their parental tumor [
1,
2]. In this respect, exploring the properties of patient-derived organoids represents an attractive option to explore tumor organization, evolution, and complex behavior in controlled conditions.
Mathematical models are often used to evaluate drug responses in infectious and oncologic disorders to predict the outcome of proposed therapy [
3,
4]. They are also used to assess the efficacy of vaccines. In this regard, nearly any signaling or metabolic pathway is describable as a combination of mathematical equations in the system [
5,
6,
7]. Previously, several authors already proposed mathematical models of tumor growth (rev. in [
8]). In such mathematical models, an ordinary differential equation (ODE) is an equation for a function of an independent variable that involves several derivatives characterizing its behavior. The systems of differential equations make it possible to address more complex problems. For instance, it would be possible to assess tumor size using the data on its cellular composition. In perspective, advanced mathematical models will define a time window most suitable for the surgery.
Moreover, several known mathematical models predict the malignancy of lung cancer. Some of these models use imaging data and the results of clinical tests (rev. in [
9]), the others – the patients' genetics, such as mutations and specific biomarkers. As believed, efficient mathematical models with good prediction performance will improve the detection of lung cancer at early stages due to the facilitation of its diagnosis.
In the present paper, we apply a computational approach to assess the parameters of a mathematical model simulating the time-dependent counts of cellular subpopulations in patient-derived organoids. The proposed mathematical model originates from the models of lung cancer progression previously developed by Geng et al. [
10]. This model accounts for four subpopulations of cancer and tumor-associated cells. Using this approach, we systematically explored changes in the cell composition of PDTOs. Here, our aim was to develop a mathematical model for prediction of changes in the cellular composition of organoids with a perspective of translating these results into the respective parental tumor which potentially can be used as a part of patient condition assessment and NSCLC prognosis.
2. Materials and Methods
2.1. Patients
The material was obtained from 16 patients from the first oncological hospital (Moscow, Russian Federation) and Novgorod Regional hospital (Veliky Novgorod, Russian federation) after surgery. The mean age was 65 years (SD 4.1), 13 of the patients were men (81%) and 3 were women (19%). All the patients have undergone thoracotomic lobectomy due to non-small cell lung cancer (NSCLC). A mathematical model was developed based on data obtained from patient-derived organoid culture of NSCLC.
2.2. Generation of PDTOs
Tumor sample was placed into basal DMEM/F12 (PanEco, Russia) culture medium supplemented with 1x penicillin-streptomycin (PanEco) and delivered to the lab on ice immediately after the surgery. In the lab, the tissue was finely chopped into small pieces using surgical scissors and then immersed in 10 mL of DMEM/F12 growth medium. This medium was supplemented with 1% Penicillin-Streptomycin (PenStrep) and 1 mg mL-1 of collagenase I (Gibco). Afterward, the minced tissue was allowed to incubate at 37°C for 1.5 hours with gentle and slow agitation.
Following the incubation, the digested tissue suspension was filtered through a 70-μm cell strainer (Falcon). The resulting cell suspension was then subjected to centrifugation at 1,500 rpm for 5 minutes at a temperature of 18–20°C. The pellet obtained after centrifugation was washed with HBSS and subsequently resuspended in 6 mL of fresh growth medium suitable for organoid culture. This organoid culture medium consisted of DMEM/F12 supplemented with 20 ng mL-1 of basic fibroblast growth factor (bFGF, 10014-HNAE, SinoBiological, China), 50 ng mL-1 of human epidermal growth factor (EGF, ab55566, Abcam), N2 supplement (PanEco, Moscow, Russia), NeuroMax (PanEco, Moscow, Russia), 10 mM Glutamax (Gibco), 1 mM N-acetyl cysteine (Sigma), 10 mM nicotinamide (Sigma), 10 μM Y27631 (ab120129, Abcam), 15 μM HEPES (Sigma), and 1% Penicillin-Streptomycin (PenStrep). The total cell count was determined using a hemocytometer, and a portion of the collected cells was preserved by freezing for subsequent flow cytometry analysis.
To generate free-floating NSCLC organoids, 96-well plates were initially coated with a 1% w/v agarose solution in Milli-Q water, applying 50 μL to each well. This agarose coating was allowed to solidify at room temperature for approximately 20–30 minutes. Subsequently, a collagen-based gel solution containing resuspended tumor cells was prepared using the SANATO 3D culture gel kit (#FTBM0051, Phystech Biomed, Russia) and 1% (v/v) SANATO reagent (#FTBM0050, Phystech Biomed, Russia) following the manufacturer's protocol. Next, 25 μL gel domes, each containing 50,000 cells, were dispensed into individual wells, and the gel was allowed to solidify at 37°C for 20 minutes. After the gel had set, 100 μL of organoid growth medium was added to each well, and the plates were placed in an incubator. The incubation was carried out at 37°C in a humidified atmosphere with 5% CO2. The growth medium was replaced every other day throughout the entire experiment. Over a period of 14 days, the organoids were subject to daily examination, during which their size and morphology were observed using bright-field microscopy with an AxioVert.A1 microscope (Zeiss, Oberkochen, Germany).
2.3. Flow cytometry
Harvested organoids were centrifuged (1,500 rpm, 5 min, r.t.). The pellet was resuspended in basal DMEM/F12 cell culture medium (1 mL/plate). Then, equal volume of 0.1% collagenase I was added and the samples were incubated for 90 min at 37°C with shacking. After the incubation, the cells were washed and resuspended in versene (0.48 mM EDTA, pH 7.4). The resuspended cells were fixed in 4% formaldehyde for 15 min at room temperature and washed twice in excess of 1X PBS. After washing the cells were resuspended in 0.1% TRITON X-100 prepared in 1x PBS and incubated for 15 min at room temperature. Then the cells were incubated with primary antibodies for 1h, in dark, and on ice. The following primary antibodies were used in this study: rabbit Alexa Fluor® 647 monoclonal anti-human antibodies directed to mannose receptor/CD206 (ab195192, Abcam), mouse monoclonal anti-human antibodies directed to CD8α (ab33786, Abcam), rabbit recombinant anti-human antibodies directed to PD-L1 (ab205921, Abcam), rabbit recombinant anti-human antibody directed to α-smooth muscle actin/αSMA (ab124964, Abcam). After the incubation, the cells were washed twice in ice-cold Versene and resuspended in 1x PBS, Then, the samples were incubated with secondary antibodies for 30 min, in dark, and on ice. The following secondary antibodies were used: goat anti-rabbit IgG H&L conjugated with Alexa Fluor® 647, preadsorbed (ab150083, USA) and goat anti-mouse IgG H&L conjugated with Alexa Fluor® 647 (ab150115, Abcam). After the incubation, the cells were washed with Versene and subjected to flow cytometry.
2.4. Mathematical approaches used in this study
The obtained results were analyzed on PC using Python 3.7. Libraries numpy 1.24.2, maptplotlib 3.7.1 were used for data processing. The scipy 1.10.1 was used to solve the differential equations system.
4. Discussion
In this paper, we analyzed changes in four subpopulations of cells in PDTOs originated from NSCLC tumors. PD-L1-positive cells represented a prominent cancer phenotype capable of suppressing the adaptive arm of immune systems. The cells expressing αSMA were typically of mesenchymal origin. The authors of experimental studies often considered αSMA as a specific biomarker of CAFs. NSCLC tumor cells positive for CD206 and CD8 represented tumor-associated macrophages and cytotoxic T-cells.
Similar to other cell-specific biomarkers, the expression of PD-L1 also occurs in healthy cells, such as immune and epithelial cells, particularly under inflammatory conditions, such as autoimmune responses, chronic infection, and sepsis [
15]. The binding of PD-L1 to the inhibitory checkpoint molecule PD-1 located on the surface of activated immune cells activates an inhibitory signal that helps the targeted cells bypass surveillance by immune cells. In this regard, the expression of PD-L1 by cancer cells lets them evade the host immune system. Recruiting stromal and immune cells allows cancer cells to bundle up an immunosuppressive microenvironment. This tumor-specific microenvironment is suitable for the proliferation of cancer cells. It also facilitates the development of drug resistance (
Figure 1).
The expression of αSMA mainly occurs in connective tissues in mesenchymal cells, such as vascular endothelial cells and smooth muscle cells. In fibroblasts, αSMA becomes induced upon their activation following an injury. In the lung, the presence of αSMA-positive fibroblasts is evident in fibrotic tissues and tumors [
16]. As a part of a tumor, the activated fibroblasts establish direct contact with the cancer cells. Moreover, cancer cells influence their gene expression by soluble factors promoting the acidification of the surrounding milieu and establishing the hypoxic conditions in the core of the tumor. Participating in the secretion and deposition of extracellular matrix [
17] by tumors, CAFs provide them with structural support and improve their resistance to anti-cancer drugs (
Figure 1).
The mannose receptor CD206 is a surface receptor of M2 macrophages. Unlike CD86-positive (M1) macrophages participating in the development of inflammatory response, M2 macrophages suppress inflammation and contribute to wound healing [
18]. The polarization of macrophages to either M1 or M2 phenotype is reversible. For instance, recruiting M1 macrophages by tumor cells induces their repolarization to the M2 phenotype [
19]. After repolarization, M2 macrophages become tumor-associated macrophages (TAM). Sustaining in the tumor microenvironment, TAMs promote tumor growth and cancer progression (
Figure 1) by producing various growth factors and cytokines, such as IL6, FGF, and VEGFA. These cytokines stimulate angiogenesis, tumor cell proliferation, and invasion. Expressing immunosuppressive cytokines (e.g., IL10 and TGFβ), TAMs suppress the activation of T-cells. TAMs also contribute to the production of extracellular matrix. In turn, covering the tumor with a layer of extracellular protects tumor cells from apoptosis and improves their resistance to anti-cancer therapies.
The CD8 protein is a specific receptor of cytotoxic T lymphocytes. This receptor participates in their interaction with antigen-presenting cells. In the body, cytotoxic CD8-positive T cells kill pathogens and eliminate cancer cells [
20]. However, some CD8-positive T-cells acquire unresponsiveness to cancer cells due to prolonged exposure to TCR signaling by cognate antigen. Chronic exposure to the antigen in the tumor microenvironment induces the genes of inhibitory receptors, such as PD-1. Although these cells remain in the tumor microenvironment, they stop expressing the proinflammatory cytokines (IFNγ, TNF, and IL2) and rarely proliferate [
21]. The exhausted CD8-positive cells contribute to tumor survival by driving tumor cells to impair the immune attack and recruiting other cells to reprogram the immune milieu. Interaction with CAFs also suppresses T lymphocyte activity [
22] by inhibiting the differentiation of CD8-positive T cells and disabling their tumor reactivity to benefit the cancer cells (
Figure 1).
Recently, other authors proposed several mathematical models linking the cellular composition of tumors with various aspects of tumorigenesis. Dean E.A. et al. [
23] constructed a mathematical model that linked the level of circulating tumor DNA with the metabolic activity of B lymphoma targeted with a T-cell-based therapy against chimeric antigen receptor (CAR). Analyzing the results of PET scans, the authors found a strong correlation between the named parameters (
r2 = 0.61,
p = 0.0001). However, they did not see an evident correlation after the exclusion of patients in remission because the metabolic activity of the tumor might reach either the maximum or upper detection limit. They also saw an even stronger correlation (
r2 = 0.79,
p = 0.0007) among the same individuals two months later due to a decline in tumor metabolic activity caused by a therapeutic response. Based on the obtained results, Dean E.A. et al. concluded that their assessments of tumor metabolic rate could be more precise if they considered heterogeneity of patients’ ctDNA.
Pedersen R.K. et al. developed a mathematical model that used measurable changes in subpopulations of tumor cells to describe the establishment of metastatic niches [
24]. Analyzing the experimental data obtained by others, the authors explored the clinical value of various processes, such as the proliferation of cancer stem cells, their differentiation, attachment, and detachment from the niche. Combining mathematical analysis and computer simulations, Pedersen R.K. et al. concluded that routine medical tests, such as blood and bulk marrow sampling, provide limited information about cellular interactions in the niche and clonal composition of participating stem cells. The authors also proposed that future studies could improve our understanding of how the metastatic niches niche form and develop and help to optimize the effects of low-dose chemotherapy.
Meade W. et al. reported a mechanistic mathematical model that predicts the effectiveness point of drugs used for prostate cancer [
25]. Analyzing the dynamics of patients’ PSA and androgen levels, the authors discovered two indicators of treatment failure. The first is the amount of androgen receptors on the surface of cancer cells, and the second is the ratio of androgen to PSA in patients’ blood. In the first case, the number of androgen receptors had to exceed a threshold value, whereas, in the second case, the ratio of androgen to PSA should fall. Both indicators predicted the loss of the therapeutic effect with high accuracy (87.3 and 88.7%, respectively). The authors concluded that using their model in clinical practice could timely advise physicians on changing the therapy.
Ait Oumghar I. et al. presented a mechanobiological bone remodeling model that couples cellular activities, mechanical stimuli, and the changes in bone mineral density during breast cancer therapy [
26]. The authors programmed and implemented the model algorithm on MATLAB software and simulated different scenarios for treatment with tamoxifen and aromatase inhibitors. Ait Oumghar I. et al. assessed their effects on bone remodeling. They also predicted how the bone volume fraction and the associated bone density loss would change over a period of time. The authors also estimated the damaging effect of each treatment regimen. Based on the results of their study, Ait Oumghar I. et al. identified the combination of chemotherapy, tamoxifen, and aromatase inhibitors, followed by the combination of chemotherapy and tamoxifen as the most risky regimen for bone damage. In addition, they confirmed this conclusion by the clinical observations of patients with breast cancer. In this regard, Ait Oumghar’s model is suitable for choosing an appropriate combination of treatments for individual patients and testing new therapeutic agents.
Larsson et al. constructed quantitative models of the time-dependent transcriptional variation of patient-derived glioblastoma cells [
27]. To build the models, the authors followed the barcoded glioblastoma cells and their progeny over three weeks. They recorded the changes in morphological characteristics and the growth rate of cancer cells. Using computational simulation, the authors discovered a hierarchical organization of the cell. They also showed that cells experienced significant phenotyping changes during their lifetime, which were somehow patient-specific. In turn, the tested therapeutic agents inhibited the specific states of the cells and altered their differentiation. In this regard, the model proposed by Larssen et al. could predict how a drug therapy affected the cancer cells in a time-dependent manner and serve for evaluation of new experimental drugs.
In turn, the mathematic model reported in this study allowed us to estimate the scale of changes in the four most abundant subpopulations of cells previously discovered in NSCLC tumors (
Table 4) and create a computer program to predict changes in cellular composition of PDTOs using the results of flow cytometry experiments performed at two data points (
Figure 2). We found that prevalence of some cells in PDTO over the others suppresses their growth (
Figure 2c and d). Cytotoxic CD8-positive cells may interfere with PDL1-positive, mostly the cancer cells (
Figure 2a and b). In
Figure 2a, the ratio of PDL1-positive cells in the sample kept declining starting from day 11 of the experiment. On the other hand, the fraction of CD8-positive cells rapidly increased. The opposite scenario occurred in PDTOs donated by patient 9 (
Figure 2b), where cytotoxic CD8-positive cells became nearly undetectable on day 14. Contrarily, the ratio of PDL1-positive cells stabilized at 3% after a steady decline. In turn, the deviations between some predicted values and experimental data indicated that the developed computer program needs improvement.
According to the previously published data, some antigens considered in our study, such as α-SMA and PDL1 (
Table 2), are not strictly limited to CAFs and cancer cells, respectively. In this regard, some cancer cells do not express PD-L1. However, expressing PD-L1 helps cancer cells to avoid immune surveillance [
28,
29]. Moreover, a certain percentage of M2 also expresses PD-L1 [
30]. Therefore, the number of PD-L1-expressing macrophages shall be subtracted from the total number of PD-L1-positive cells. In turn, some cancer cells express α-SMA because α-SMA is one of the known biomarkers of epithelial-mesenchymal transition. However, the number of non-fibroblast α-SMA-positive cells in our experimental conditions was negligible by measuring the cellular composition in the experimental block. In addition to the above, a systematic literature review was conducted throughout the model study pathway to more accurately reproduce intercellular interactions. Quantification of double-labeled cells and taking them into account would improve the accuracy of the made predictions.
In turn, the presented mathematic model also has several limitations. First, we performed this study on a limited number of patients. Respectively, a higher sample size would improve the accuracy of our predictions. Second, our prediction power is time-limited since we used the experimental data obtained at three different time points - on days 7, and 14 of the experiment. Respectively, getting out of this time interval will lower the accuracy of our assessments. Third, although organoids are similar to their parental tumors in many aspects [
31], their prolonged culturing in vitro would result in different clinical phenotypes due to the isolation of organoids from the host immune and endocrine system. The organoids do not experience pressure from the approaching immune cells. The growth factors that the organoids receive with culture medium are not necessarily the same and in the same concentrations as the ones delivered to the tumor in vivo. [
32]. Fourth, the ability of organoids to recruit new cells is limited to the cells already present in the well at plating. Fifth, the model does not consider changes in minor subpopulations of cells.