An agent-based model of monocyte differentiation into tumour-associated macrophages in chronic lymphocytic leukemia

Abstract

Monocyte-derived macrophages help maintain tissue homeostasis and defend the organism against pathogens. In tumors, recent studies have uncovered complex macrophage populations, including tumor-associated macrophages, which support tumorigenesis through cancer hallmarks such as immunosuppression, angiogenesis, or matrix remodeling. In the case of chronic lymphocytic leukemia, these macrophages are known as nurse-like cells (NLCs) and they protect leukemic cells from spontaneous apoptosis, contributing to their chemoresistance. We propose an agent-based model of monocyte differentiation into NLCs upon contact with leukemic B cells in vitro. We performed patient-specific model optimization using cultures of peripheral blood mononuclear cells from patients. Using our model, we were able to reproduce the temporal survival dynamics of cancer cells in a patient-specific manner and to identify patient groups related to distinct macrophage phenotypes. Our results show a potentially important role of phagocytosis in the polarization process of NLCs and in promoting cancer cells' enhanced survival.

Keywords: Cancer; Computational bioinformatics; Health sciences; Immunology.

Graphical abstract

Figure 1

Experimental setups and datasets from in vitro PBMC cultures from CLL patients (A) Experimental set-up. Autologous and heterologous co-cultures of CLL cells and monocytes leading to NLC formation. In autologous cultures, Peripheral Blood Mononuclear Cells (PBMC) were isolated from CLL patients’ blood samples and cultured in vitro for 13 days. The cell concentration and viability of CLL cells was monitored by hemocytometer and flow cytometry AnnexinV/7AAD staining, respectively. (B) Visualization of NLCs at 10 days of in vitro culture from two different patients in bright-field and immunofluorescence microscopy (NLC: green staining; CLL cells: red staining). (C) Time course datasets produced from the PBMC autologous cultures from 9 patients. CLL cell survival was monitored by viability assay and concentration measurements. The black curve corresponds to the mean value averaged over the available data. The shaded area corresponds to the 95% confidence interval. The complete dataset showing patients variability is available in Supplementary Material (Figure S2). Time points at which data was not available for at least 4 patients were removed from downstream analysis (Day 4, 5, 11, and 12). (D) Heterologous co-cultures. Monocytes from healthy donors and B cells from CLL patients were co-cultured to assess the relationship between the initial density of monocytes in the culture and the level of survival of CLL cells after 9 days. The x axis displays different monocytes initial proportions (not to scale for clarity). Measurements were performed on co-cultures of B cells from 5 CLL patients and monocytes from 2 healthy donors. Boxplots show the interquartile interval (25–75%) and the median of the measurements for each initial density of monocytes. Whiskers extend to the furthest datapoint within the 1.5x interquartile range and more extreme points are marked as outliers. The complete data showing inter-patient and inter-donor variability is available in Supplementary Material (Figure S3).

Figure 2

ABM representations (A) Netlogo simulation of 5000 cells. Cancer cells are depicted as small arrows (red, yellow or gray for NeedSignal, Apoptotic or Dead state, respectively) and myeloid cells are depicted as pentagons (blue, orange or green for Monocyte, Macrophage or NLC, respectively). (B) Schematic diagram of the agents’ states and behaviors. Parameters optimized through the genetic algorithm are indicated next to their corresponding cellular processes, represented by arrows.

Figure 3

Overview of strategy and results of the parameter exploration (A) Schematic diagram of the biobjective genetic algorithm. Flow chart of the NSGA-II algorithm procedure. The algorithm parameters and procedure details are described in the STAR Methods section. (B) Pareto front of the biobjective optimization. The NSGA-II genetic algorithm evaluates each explored set of parameters according to 2 objective functions corresponding to CLL cell viability and concentration (10 time points, least squares method), along 20′000 generations. The Pareto front, depicted in red, contains 98 non-dominated solutions. The knee-point parameter set is in a cyan box and the parameter sets performing best on viability and concentration dynamics are in magenta boxes. (C) Representative parameter sets. The parameters listed here correspond to the knee-point parameter set, and to the parameter sets fitting best viability and concentration dynamics.

Figure 4

Parameter distributions and correlations (A) Parameter distribution of the searched space and of the Pareto front solutions. The violin plots show the distribution of parameter values in the entire search space after parameter exploration throughout 20′000 generations. The swarm plots represent the 98 non-dominated solutions of the Pareto front. The cyan dot in each plot corresponds to the corresponding parameter value in the knee-point set. (B) Parameter correlations with fitness on viability. The parameters are ranked based on the absolute value of their correlation coefficient with fitness on viability. (C) Parameter correlations with fitness on concentration. The parameters are ranked based on the absolute value of their correlation coefficient with fitness on concentration.

Figure 5

Comparison of simulated and experimental results (A–C) Model fitting on PBMC autologous cultures. 12 simulations were run with the knee-point parameter set (A), the parameter set maximizing the viability fitness (B), the parameter set maximizing the concentration fitness (C), and compared with the experimentally observed viability and concentration dynamics averaged over 9 patients. The initial monocyte and apoptotic cancer cell proportions for the simulations were set to the average monocyte and apoptotic cell proportions measured in the patient PBMCs (1.28% and 4.32%, respectively). Simulations are depicted in red and experimental data in black. The black curve corresponds to the mean value averaged over the available data. The shaded area corresponds to the 95% confidence interval. (D) Model predictions on heterologous co-cultures with varying monocyte initial proportions. Simulations were run varying initial monocyte proportions (3 repetitions) for 9 days and are here compared to experimental observations in heterologous co-cultures with the corresponding initial conditions after 9 days (average over 10 experiments including 5 CLL patients and 2 healthy donors). Red dots correspond to the simulations and black dots and boxes correspond to the experimental data. Boxplots show the interquartile interval (25–75%) and the median of the measurements for each initial density of monocytes. Whiskers extend to the furthest datapoint within the 1.5x interquartile range and more extreme points are marked as outliers. Values of ${\mathrm{R}}^2$ approaching one and values of NRMSE close to zero indicate a good performance of the model.

Figure 6

Parameter sensitivity analysis The parameters were varied one at a time while keeping all other parameters constant to estimate their impact on the overall dynamics (4 simulation runs per value). Parameters which had the largest impact are displayed here (based on having absolute correlation coefficient to fitness on viability or on concentration above 0.4). A parameter sensitivity analysis on the remaining 13 parameters is shown in Figure S8. The black curve corresponds to the mean value averaged over the available experimental data. The shaded area corresponds to the 95% confidence interval.

Figure 7

Evaluation of the patient-specific models (A) Model fitting of patient-specific models on PBMC autologous cultures. On each panel, 12 simulations are shown with the corresponding patient-specific knee-point parameter set and compared with the experimentally observed viability and concentration dynamics. The initial monocyte proportion for the simulations was set to the corresponding monocyte proportion measured in each patient (Table S2). Simulations are depicted in red and experimental data in black. (B) Prediction performances of 3 patient-specific models on heterologous co-cultures with varying monocyte initial proportions. Simulations were run for varying initial monocyte proportions (3 repetitions) for 9 days and are here compared to experimental observations in heterologous co-cultures with the corresponding initial conditions after 9 days. Red dots correspond to the simulations and black dots correspond to the experimental data. Values of ${\mathrm{R}}^2$ approaching one and values of NRMSE close to zero indicate a good performance of the model. For each patient, experiments were carried out with varying proportions of monocytes from 2 different healthy donors. However, due to low sample quantities from either the patient and/or the donor, not all monocyte proportions could be tested for all patients. The complete data showing inter-patient and inter-donor variability is available in Supplementary Material (Figure S3).

Figure 8

Patient clusters (A) Unsupervised hierarchical patient clustering. Patients were clustered with the complete linkage clustering method based on their knee-point parameters sets. Parameter values were centered and scaled by row. Patient-specific parameter sets are available in detail in Table S3. (B) Principal Component Analysis. Patients are colored according to the class identified in the unsupervised hierarchical clustering shown in A. (C) Parameter contributions to the first dimension of the principal component. The dashed line corresponds to the expected value if the contributions were uniform. Any parameter with a contribution above the reference line could be considered as important in contributing to the dimension.

Figure 9

Comparative knee-point parameter sets distributions within each patient class The knee-point parameter sets of patients from clusters A and B resulting from the unsupervised clustering shown on Figure 8A were integrated and depicted here in blue for Class A and orange for Class B. Boxplots show the interquartile interval (25–75%) and the median of the parameters’ values. Whiskers extend to the furthest datapoint within the 1.5x interquartile range and more extreme points are marked as outliers.

Figure 10

Unsupervised patient classification reveals characteristic macrophage profiles. The black curve corresponds to the mean value averaged over the available data for each class. The shaded area corresponds to the 95% confidence interval.