A phase-field model for non-small cell lung cancer under the effects of immunotherapy

1.   Introduction

A major goal of mathematical oncology is predicting the growth of tumors ^[1]. Cancer is a class of diseases characterized by numerous point mutations in the genome that result in the uncontrolled growth and spread of cells. Overviews of its biological details are, e.g., presented in ^[2,3]. The body's immune system can suppress tumor growth by inhibiting cell growth or by destroying cancer cells. On the other hand, it can also promote tumor progression by selecting tumor cells that are better able to survive in an immunocompetent host or by establishing conditions within the tumor microenvironment that facilitate tumor outgrowth ^[4,5]. Immunotherapy attempts to boost the body's immune system and its immune responses against cancerous cells to eliminate them. Understanding this ability has revolutionized the field of oncology ^[6]. For a variety of cancer types, immunotherapy has been proven to provide a significant clinical benefit in patients with advanced stages of cancer and is, as of today, well-established as a standard treatment method ^[7,8,9]. However, it is extremely difficult to accurately identify spatial tumor growth on a patient-specific level, especially under a treatment plan ^[10,11]. A variety of mathematical models has helped to improve the understanding of biological principles in cancer growth and treatment response ^{[12,13,14,15,16]} and their predictive power ^[17,18]. Globally, lung cancer is the leading cancer-related cause of mortality. Roughly 85 percent of all lung cancer cases are non-small cell lung cancer (NSCLC). Its five year survival probability is approximately 22 percent ^[19]. Current treatment options for NSCLC include surgery, but also chemotherapy, radiotherapy, and, in recent years, immunotherapy ^[20]. NSCLC can be again divided into three different major subgroups: Adenocarcinoma, squamous cell carcinoma, and large cell carcinoma. Immunohistological testing can identify the respective subgroups to which a patient's tumor belongs ^[21].

For overviews of spatial tumor growth models, including continuum, discrete, and hybrid models, we refer to the reviews in ^[22,23]. One class of continuum models relies on simple reaction-diffusion-equations and elastic models for mass effects ^[24,25,26]. We will rely on phase-field models commonly used for modeling cell dynamics, since they allow us to model the observed cell-to-cell adhesion ^[27] between tumor cells by energy terms. Simple two-species models consisting of a phase-field model for the tumor and a nutrient equation of reaction-diffusion type are introduced, analyzed or verified in ^[28,29,30]. Models, including increasingly more complicated flow models, are given in ^{[31,32,33,34,35]}. A more general theoretical approach to multispecies models can be found in ^{[36,37,38,39,40]}. The number of modeled species and fields strongly depends on the choice of the particular studied effects. More specialized large models can be found in ^{[41,42,43,44,45]}.

Ordinary-differential equation-based models used to simulate the immune response of cancer cells, as well as simple immunotherapy models, are described in ^{[46,47,48,49]}. Spatial models try to account for a variation of the tumor structure in the entire spatial domain. Regarding those models, which also include cancer therapy, Powathil et al. ^[50] introduces a hybrid model to study the impact of different chemo- and radiation therapies on tumor cells. The chemotherapy of breast cancer with a drug-delivery model is discussed in Wu et al. ^[51], and with an additional reaction-diffusion model of tumor growth in ^[35,52]. In Hormuth et al. ^[53], a continuum model involving radiation therapy for brain tumors is discussed. For prostate cancer, a combination of chemotherapy with antiangiogenic therapy and the optimization of the treatment protocols is given in Colli et al. ^[54]. Chemotherapy is also included in the multispecies phase-field model approach of Fritz et al. ^[55].

The present work aims to develop a mathematical model for the spatial growth behavior of solid tumors in NSCLC in a simple, deterministic setting. It addresses the effects of immunotherapy application and allows the description of its influences on the tumor's spatial structure. The model framework is applied to two data sets acquired from clinical patients that have shown qualitatively different therapy outcomes. Analysis of the model sheds light on the corresponding parameter relations that determine different clinical outcomes that could potentially be estimated in a prognostic setting to improve the prediction of clinical outcomes and therapy choices.

Here, passive immunotherapy that uses monoclonal antibodies to improve the immune response by regulating T-cell activity in the effector phase of the immune response via the so-called PD-1/PD-L1-pathway is considered. The downregulation usually caused by PD-1 activation prevents collateral damage during an immune response ^[56,57]. However, tumor cells can interfere with this pathway by expressing the corresponding ligands PD-L1 and PD-L2 that bind to the T-cells PD-1 receptor and inactivate the T-cell to decrease the immune response towards the tumor cells ^{[56,58,59,60]}. Cells expressing these ligands are targets for, e.g., the drugs Nivolumab and Pembrolizumab with which patients of our clinical data set were treated.

2.   Models

We consider a bounded domain $\Omega \subset \mathbb{R}^3$ representing a volume of lung tissue in which a cluster of tumor cells reside, their volume fraction given by a function $\phi_T$ of position $x \in \Omega$ and time $t > 0$. We divide $\phi_T$ into two types of cell species, of proliferative tumor cells $\phi_P$, and necrotic tumor cells $\phi_N$, such that $\phi_T = \phi_P + \phi_N$. We do not explicitly model hypoxic cell species or differentiate between immunotherapy sensitive or resistant cells as in the references ^[49,61]. Note, that by introducing the healthy cells $\phi_H$, we have the partition $1 = \phi_P + \phi_H + \phi_N$ of volume fractions. Furthermore, we introduce the nutrient concentrations $\phi_{\sigma, v}$ and $\phi_{\sigma, i}$ inside the vasculature and interstitial space. The concentration of immunotherapeutic agents is given by $\phi_\tau$. We will first start with the full 3D-Model and then introduce its simplifications to a 1D-Model.

2.1.   The full 3D-Model

We use a generalized Cahn–Hilliard model as in ^[41,44] with reaction terms to describe the evolution of proliferative cells. By introducing the chemical potential $\mu_P$, the model is characterized on the computational domain $\Omega$ by the fourth-order equation system

with boundary conditions

Here, $\Psi(\phi_P, \phi_T) = \tilde{\Psi}(\phi_P) + \tilde{\Psi}(\phi_T)$, with $\tilde{\Psi}(\phi) = c_\Psi \phi^2 (1-\phi)^2$ is a double-well potential with scaling factor $c_\Psi$ that incentivises $\phi_P$ and $\phi_T$ to become either 0 or 1. The interfacial width of our phase-field is given by $\varepsilon_P > 0$, which scales the kinetic energy term. The constant $\chi > 0$ describes the strength of the chemotaxis induced by the nutrients on the tumor cells ^[62]. Finally, the cell mobility of the tumor cells is given by $m_P(\phi_P, \phi_T) = (\phi_P^{+})^2 (1-\phi_T^{+})^2$ with the constant $c_m > 0$. Here, $(\cdot)^+$ denotes a cutoff function restricting fields to values between zero and one. The mobility is designed to force the tumor cells to stay between 0 and 1. We remark that the boundary conditions are chosen to support mass conservation in $\Omega$ in the absence of the reaction term. The reaction term, $S_P$, is given by

The first term on the right-hand side of (2.2), is a Gompertzian growth term ^[63] acting on the tumor interface, with a small parameter $\epsilon_{g}$ for regularization. We assume that the growth is proportional to the nutrient concentration and scale it by the parameter $\lambda^{pro}_P$. For increasing $\lambda \geq 1$, the growth of smaller tumor cell concentrations is penalized, which inhibits the tumor from spreading over the entire lung.

The second term on the right-hand side of (2.2) models the decay of tumor cells due to immunotherapy and acts on the tumor volume. Here, $\frac{1}{\omega}$ is the ratio of the lung weight to the body weight, $\lambda^{ \rm{eff}}_\tau$ describes the patient-specific effect of the immunotherapy on the tumor, and $\phi^{50}_\tau$ is the drug concentration for the half-maximal response.

The last term of (2.2) models a decay of proliferative cells due to low nutrient concentrations with rate $\lambda_{PN}$. It is activated by the Heaviside function, $\mathcal{H}$, if the nutrient concentration becomes smaller than $\sigma_{PN}$.

The necrotic cell field is described by a simple spatial ODE model:

where the growth term on the right-hand side mirrors the decay terms of the proliferative cells with a different sign. In the absence of growth ($\lambda_P^{pro} = 0$) and therapy ($\lambda^{ \rm{eff}}_\tau = 0$) and, due to the boundary conditions mass conservation for the entire tumor $\phi_T$ is obtained.

For the nutrients, we follow the approach in Erbertseder et al. ^[64] and model their distribution with a stationary diffusion model:

and

with homogeneous Neumann boundary conditions: $\frac{\partial \phi_{\sigma, v}}{\partial n} = 0$ and $\frac{\partial \phi_{\sigma, i}}{\partial n} = 0$ on $\partial \Omega$. Static nutrient models are used in ^[65,66] and are based on the observation that the tumor evolves on a timescale of days and weeks, while the nutrient model works on a much faster time scale. Here $- \xi_{va} \phi_{\sigma, v}$ models the coupling with the arteries, while $- \eta_{vi} \phi_{\sigma, v}$ and $\eta_{iv} \phi_{\sigma, i}$ model the nutrient exchange between the vasculature and interstitial space. The source term couples the equation to a 2D manifold $\Gamma$, approximating the surface of arteries and bronchial tubes, by a Delta-Distribution. In our case, the geometric information of $\Gamma$ originates from patient-specific CT scans from the beginning of the treatment period. The amount of nutrients entering the system from outside is controlled by the parameter $\beta > 0$. The nutrient consumption by healthy cells $\phi_H = 1 - \phi_P - \phi_N$ is modeled by the decay term $-\alpha_H (1 - \phi_P - \phi_N) \phi_{\sigma, i}$, while $-\alpha_P \phi_P \phi_{\sigma, i}$ models the enhanced consumption by proliferative cells. Finally, $\kappa_v$ and $\kappa_i$ are the vessel's and interstitial space's scalar diffusion constants.

We model the immunotherapy concentration with the ODE,

where the first term describes the decay with the drug serum half-life time $t_{1/2}$, and the second term is the drug influx due to injections. The treatment plan is described by the set $I_\tau \subset [0, \infty)$, which contains the times at which the drug is administered. Hence, the indicator function acts as an activation function for the source term. Furthermore, $N_A = 6.022140857\cdot 10^{23}$ denotes the Avogadro constant, $M_\tau$ is the molar mass of the drug, and $d_\tau = 0.24$ is the administered dosage of the antibodies. As done in ^[16], we assume that the immunotherapeutic concentration follows the behavior of a Hill–Langmuir equation. All the model parameters are summarized on the left of Table 1.

Table 1.  List of model parameters. Left columns: Mathematical symbol. Center columns: Dimension and description. Right columns: Numerical value for our simulation and references.

	Dim.	Description	Values 1D		3D
			Patient 1	Patient 2	Patient 1
$c_{m}$	$\left[m^2/d \right]$	Mobility factor of proliferative tumor cells.	$5\cdot10^{-4}$			[67]
$\lambda_P^{\rm{pro}}$	$\left[1/d \right]$	Proliferation rate of proliferative tumor cells.	0.36	0.0038	0.47
$\lambda$	$\left[\right]$	Power-law for growth of cell concentrations.	2	1	2
$\lambda_{PN}$	$\left[1/d \right]$	Decay rate proliferative to necrotic cells.	0.1	20	0.05
$\sigma_{PN}$	$\left[\right]$	Nutrient threshold for necrotic growth.	0.2
$\varepsilon_P$	$\left[m \right]$	Tumor interface width.	$5\cdot10^{-4}$
$c_\Psi$	$\left[\right]$	Scaling double-well potential.	2
$\epsilon_g$	$\left[\right]$	Regularization for Gompertzian growth term.	0.1
$\chi$	$\left[\right]$	Chemotaxis factor.	0
$\kappa_v$	$\left[m^2/d \right]$	Nutrient diffusion factor in the vasculature.	-		$10^{-3}$
$\kappa_i$	$\left[m^2/d \right]$	Nutrient diffusion factor in the interstitial space.	$10^{-5}$
$\xi_{va}$	$\left[1/d \right]$	Coupling smaller to larger arteries.	-		0
$\eta_{vi}$	$\left[1/d \right]$	Coupling vasculature to interstitial space.	1
$\eta_{iv}$	$\left[1/d \right]$	Coupling interstitial space to vasculature.	0
$\alpha_H$	$\left[1/d \right]$	Nutrient consumption rate healthy cells.	1
$\alpha_P$	$\left[1/d \right]$	Nutrient consumption rate proliferative cells.	4
$\beta$	$\left[m/d \right]$	Source term large vessels.	-		0.1
$\lambda_\tau^\text{eff}$	$\left[1/d \right]$	Immunotherapeutic effect under application.	5.2	0.54	4.49
$\phi_\tau^{50}$	$\left[mol \right]$	Drug concentration for half-maximal response.	$1.012\cdot 10^{16}$			[16]
$d_\tau$	$\left[g \right]$	Dosage per medication interval.	0.24	0.20	0.24	[16]
$M_\tau$	$\left[g/mol \right]$	Molar mass of medication.	146.000	143.600	146.000	[16]
$t_{1/2}$	$\left[d \right]$	Drug serum half-life time.	26.7	22.0	26.7	[16]
$\omega$	$\left[\right]$	Ratio of the patient weight to the lung weight.	80

---

 | Show Table

DownLoad: CSV

2.2.   The reduced 1D-Model

Since the full 3D model has high computational costs, we begin by assuming a spherically symmetric tumor to precalibrate our model parameters. Thus, we begin with a model in which we assume a radial dependency of $\phi_P = \phi_P(r)$ and $\phi_{\sigma, v} = \phi_{\sigma, v}(r)$, which immediately implies $\mu_P = \mu_P(r)$, $\phi_{\sigma, i} = \phi_{\sigma, i}(r)$ and $\phi_N = \phi_N(r)$. As a first step, we assume spherical symmetry and, hence, obtain the following simplified model in the radial coordinates:

with boundary conditions

The equations for the necrotic cells and immunotherapeutic agents stay unchanged. For the nutrient model, we assume that the nutrient concentration $\phi_{\sigma, v}$ can be approximated by a pre-given constant value and use

for the nutrients in the interstitial space.

We remark that the $r^{-2}$ term does not pose any numerical issues for finite elements since we have to apply the spherical volume measure $4 \pi r^2 \text{d}r$ when we bring the system into its weak form and, thus, the $r$-terms cancel.

3.   Methods & Materials

Our numerical discretization in space and time consists of finite elements and a semi-implicit scheme in time. The lung geometry in our model is taken from CT data obtained from two patients. Their therapy plans are presented, and the measured tumor data is introduced. Further details follow.

3.1.   Numerical discretization

As a spatial discretization, we use finite elements with piecewise-linear basis functions for $\phi_P$, $\mu_P$, $\phi_N$, $\phi_{\sigma, v}$, and $\phi_{\sigma, i}$. We use a semi-implicit scheme in time for (2.1). The source and sink terms, as well as the mobilities, are treated explicitly. For the potential $\Psi$, we use a simple convex-concave splitting in time: Note that the double-well potential $\tilde\Psi$ can be decomposed into the $\tilde\Psi_i = \frac{3}{2} c_\Psi \phi^2$ and $\tilde\Psi_e = c_\Psi \left(- 2\phi^3 -\frac{1}{2} \phi^2 + \phi^4 \right)$. Clearly, $\Psi_i$ is convex, while $\Psi_e$ stays concave for $\phi \in \left[\frac{1}{2}-\frac{1}{\sqrt{3}}, \frac{1}{2}+\frac{1}{\sqrt{3}} \right] \supset \left[0, 1\right]$. In our case, this motivates the decomposition $\Psi_i (\phi_P, \phi_N) = \tilde{\Psi}_i(\phi_P) + \tilde{\Psi}_i(\phi_P + \phi_N)$ and $\Psi_e (\phi_P, \phi_N) = \tilde{\Psi}_e(\phi_P) + \tilde{\Psi}_e(\phi_P + \phi_N)$. For a fixed $\phi_N$, the latter is concave if $\phi_P \in \left[\frac{1}{2}-\frac{1}{\sqrt{3}}, \frac{1}{2}+\frac{1}{\sqrt{3}} - \phi_N \right] \supset \left[0, 1- \phi_N \right]$. In the case of no source terms, we achieve unconditional gradient stability ^[68] if $\Psi_i$ is treated implicitly and $\Psi_e$ explicitly and the bounds of $\phi_P$ are satisfied. As a linear solver, we use MINRES with the block-diagonal preconditioner

which is a generalization to the preconditioner in ^[69]. For the inversion of the diagonal blocks, we resort to algebraic multigrid ^[70].

For (2.3) we use an explicit Euler scheme, while in (2.4) and (2.5) $\phi_P$ is kept explicit, and we solve for $\phi_{\sigma, v}$ and $\phi_{\sigma, i}$. Since in our studies, we set $\eta_{iv} = 0$, both equations are decoupled and can be solved separately with a Conjugate Gradient method preconditioned by algebraic multigrid ^[70].

Finally, we apply a simple explicit Euler scheme for the ODE model of (2.6). Since (2.6) runs physically on a much smaller time scale, we implement it with a time step which is by a factor of ${}^{1}\!\!\diagup\!\!{}_{n}\;$ smaller than the time step of (2.3). The coupling does not pose difficulties since (2.6) is completely decoupled and computationally cheap to solve. For our implementation, we use the FEniCS framework ^[71].

3.2.   Simulation geometry

The lung geometry and vasculature structure were extracted from the patient's CT data (see Patient Data) at the therapy's onset. For the lung geometry, the extraction was performed with 3D Slicer ^[72], which was simplified with 3D Builder by Microsoft Corporation and blender ^[73] before remeshing it coarsely with gmsh ^[74]. The vasculature was extracted using slicers' vmtk extension ^[75,76].

The vasculature and the lung are depicted in Figure 1(A). To keep the computational costs tractable, we apply a local refinement close to the initial tumor position by subdividing the mesh several times, see Figure 1(B). This approach is motivated by the form of the potential energy of the Cahn–Hilliard equation, which keeps the central part of the tumor localized to its initial location. Similarly, we refine the mesh close to the tumor for the nutrient model, which couples with the Cahn–Hilliard equation for improved accuracy.

Figure 1.  Patient-specific lung geometry reconstructed from CT images: (A): Lung envelope and 2D-surface of the vasculature. (B): 2D slice of the 3D mesh near the tumor.

DownLoad: Full-Size Img PowerPoint

3.3.   Data

Datasets from two clinical patients were available to support our model as a proof of concept. The examined patient data consist of anonymized volumetric measurements of primary tumors and metastases in patients with NSCLC and was acquired from routinely generated CT slices and computations of the secondary appraisal environment syngo.CT Lung Computer Aided Detection (LCAD) workflow of Siemens Healthineers, provided in the syngo.Via framework.^* As measurements, we have the manually segmented tumor volume, the response evaluation criteria in solid tumors (RECIST v1.1, ^[77]), and the bidimensional World Health Organization criteria (WHO), all indicating the tumor response with respect to therapy in different defined ways. Since the LCAD environment identifies suspicious nodules below the usual detection size threshold on CT scans as metastases, we did consider these in the measurements to provide a more physically accurate and acquisition-consistent representation.

^*The patients were treated in the Clinic of Pneumology, Thoracic Oncology, Sleep and Respiratory Critical Care of the Klinikverbund Allgäu in Germany. The ethics commission of BLAEK (Ethik-Kommission der Bayerischen Landesärztekammer), approved the use of the data (reference number 19021).

3.3.1.   Patient 1

The patient data during the therapy are depicted in Figure 2(A), which contains the tumor volume and the RECIST and WHO values. The tumor was diagnosed at time $t = 0$, and the therapy started at time $t_{S} = 296$ days. The patient received Nivolumab as an immunotherapeutic antibody every two weeks, for 30-60 minutes, in an intravenous dose administration from this time on. At time point (E), the patient developed a temporary lung emboly, significantly influencing the data quality and LCAD accuracy. From (P) on, a break of immunotherapeutic drug administration was conducted due to the evident presence of brain metastases. At that very time point, the treating physicians diagnosed a complete remission for the primary tumor in the lung. All in all, we model our therapy as

Figure 2.  Measured patient data from CT scans. Top row: tumor volume. Middle row: RECIST value. Bottom row: WHO value. Vertical dashed lines indicate events that might affect the data quality or therapy. For Patient 1 E indicates an Embolie, P a therapy break, and D the treatment with Dexamethason. For Patient 2 Q3W is the beginning of a three-weekly drug administration schedule, Q6W the change to a six-weekly administration, PD the therapy end, and EL the Exitus Letalis.

DownLoad: Full-Size Img PowerPoint

3.3.2.   Patient 2

For the second patient, the therapy starts before we have our first data point. The patient receives 200 mg of Pembrolizumab every three weeks. The therapy plan changes at time $t_{Q6W} = 514$ days, where the application interval is adjusted from three to six weeks. Finally, it was terminated at time $t_{PD} = 1063$ days, followed by the patient's death (exitus letalis) after $1306$ days. The final CT image shows a very large tumor mass that has spread over the entire lung. The corresponding therapy for this patient is modeled as

4.   Results

In this section, we discuss our simulation results. All the used parameter values are summarized on the right of Table 1. In CT images, only a small part of the tumor cells are visible. To take this into account, we define the visible tumor volume from our simulation as $V_{vis} = \int_\Omega \mathbf{1}_{\{x \in \Omega : \phi_{T}(x) > 0.3\}}(x) \text{d}x$. Numerically, we apply a time step size of ${}^{1}\!\!\diagup\!\!{}_{24}\;$ for the tumor and $n = 32$ substeps for the evolution of the immunotherapeutic agents. For the 1D model, we use spatially 500 elements in a spherical domain with radius 4 cm.

4.1.   1D results

We first discuss the results obtained with the simplified 1D model in spherical coordinates. Figure 3(A) shows a comparison of the results obtained for the first patient. The manually calibrated parameters are $\lambda^{ \rm{pro}}_P = 0.36, \lambda^{\text{eff}}_\tau = 5.2$, $\lambda = 2$ and $\lambda_{PN} = 0.1$. We start with an initial tumor volume with enough temporal distance from the previous therapy window. At the top, the visible tumor volume of our simulation is shown as a dotted line, which qualitatively follows the segmented data points. The number of immunotherapeutic agents is depicted in the middle. Note that after a few months, a quasi-periodic state is reached. The lower plot depicts the entire tumor mass in the whole domain and therefore includes cell concentrations that cannot be detected yet. Its evolution mainly depends on the parameter $\lambda$, such that larger values of $\lambda$ impede the tumor cells from spreading over the entire domain. In this case, the tumor mass curve follows the tumor volume curve. All in all, the given immunotherapy manages to eliminate the proliferative tumor cells. Since, at the end of the therapy only necrotic cells remain, the tumor stays under control, even though no immunotherapeutic agents hamper its growth.

Figure 3.  Patient data and simulation results: Contains the visible tumor volume, the number of immunotherapy agents, and the tumor mass over time. (A): Simulation results for Patient 1. (B): Simulation results for Patient 2.

DownLoad: Full-Size Img PowerPoint

Figure 4 shows the simulation state in spherical coordinates. In Figure 4(A), the initial state of the tumor simulation is shown. The entire tumor consists of proliferative cells, and no necrotic cells have been able to form yet. We see a large drop in the nutrient concentration due to the increased consumption of the tumor cells. In Figure 4(B), we observe the tumor volume fraction just before the start of the immunotherapeutic treatment. Because of the scarcity of nutrients at the tumor core, some of the inner proliferative cells have necrotized. Since necrotic cells do not deplete nutrients, the nutrient concentrations do not decrease even though the entire tumor has grown. Figure 4(C) is a snapshot taken during the therapy. Since the immunotherapeutic term acts on the entire volume of proliferative cells, while the growth terms only act on the interfaces, the proliferative cells have increased their interface to the other cell species by becoming active at the tumor boundary. This leads to an unexpected increase of nutrients near the necrotic tumor core. In Figure 4(D), the proliferative tumor shell contracts due to the therapy until the necrotic and proliferative cells form a small tumor sphere. Due to a large nutrient concentration, no conversion of proliferative to necrotic cells happens at this point. The ongoing immunotherapeutic treatment decreases the proliferative cells until, in the end, only the already present necrotic core remains, which is in agreement with our clinical observations.

Figure 4.  Cell- and nutrient fields for Patient 1 plotted over radial distance. (A): Start of the simulation. (B): Just before the therapy. (C): Just after therapy starts. (D): During the therapy.

DownLoad: Full-Size Img PowerPoint

Figure 3(B) shows a comparison of the results obtained for the second patient. Here, the therapy is already in progress, and different from the first patient, we have no time window to observe the tumor without any therapy. The manually calibrated parameters are $\lambda^{ \rm{pro}}_P = 0.0038, \lambda^{\text{eff}}_\tau = 0.54$, $\lambda = 1$ and $\lambda_{PN} = 20$. With respect to the tumor volumes in the topmost plot, the first data point is ignored and a plausible tumor volume is assumed, more in line with the previous values. The visible tumor diminishes in size until it is able to sustain itself with the given nutrients despite ongoing therapy and stays constant. The effect to the change in the injection schedule from 3 to 6 weeks shows itself in a decreased amount of immunotherapeutic agents, but has no notable impact on the visible tumor. Even though the visible tumor does not grow, its cells start to spread over the entire lung, which leads to a consistent increase in the total tumor mass in the lower plot. After the tumor has spread over the entire domain, the visible tumor starts to grow rapidly, which is in line with the observed outcome of the tumor patient.

Figure 5 gives a more detailed view of what happens spatially. Figure 5(A) describes a similar initial setup to the first patient: The entire tumor consists of proliferative cells. No necrotic core has formed yet, but the nutrient concentration is notably lower inside the tumor. After 170 days, the state illustrated in Figure 5(B) is reached, where a necrotic core has formed, the proliferative cells have developed an outer shell, and the nutrient concentration inside the tumor has stabilized. We also observe that a small amount of proliferative cells has moved away from the primary tumor and spread towards the outer boundary. In Figure 5(C), after 600 days, the proliferative shell has collapsed, and we have a tiny amount of proliferative cells which can sustain themselves due to the increased nutrient concentration. This configuration is stable over the next year, as we see from Figure 5(D), which is 1175 days after the start of the simulation. This is just before the visible tumor starts to grow again. We observe that the tumor has spread over the entire domain, which leads to a visible decrease in nutrient concentration from 1 to roughly 0.9. At this point, the phase separation from the Cahn–Hilliard is initiated, resulting in drastic tumor growth within a short period. We point out that, especially in the case of explosive growth in the second patient, the simulation geometry does play a role. This is the main motivation behind moving from the simplified model in spherical coordinates to more elaborate 3D simulations.

Figure 5.  Cell- and nutrient fields for Patient 2 plotted over radial distance. (A): Simulation start. (B): After 170 days. (C): After 600 days. (D): After 1175 days.

DownLoad: Full-Size Img PowerPoint

Under the assumption that we have meaningful parameter values for both patients, an interesting follow-up question is to investigate how changes in the therapy qualitatively affect its outcomes.

4.2.   Different therapy plans

Figure 6 shows the tumor evolution for both patients with changed therapy plans. We monitor the visible tumor volume at the top, the immunotherapeutic drug concentration in the middle, and the total tumor mass, including unobservable components at the bottom. Contributions due to necrotic cells are indicated with dashed lines, proliferative tumor cells have dotted lines, immunotherapeutic drugs are depicted with solid lines, and the visible parts of the tumor are assigned dash-dotted lines. The experimental data is always depicted in black color, while the original therapy plan is depicted in blue. For the alternative therapy plans red, green, and cyan colored curves are used.

Figure 6.  Patient data and simulation results: Contains the visible tumor volume, the number of immunotherapy agents, and the tumor mass over time for different therapies. Experimental data is depicted in black, the original therapy is depicted in blue, while alternative therapy plans are given in green, red, and cyan. (A): Different drug dosages for Patient 1. (B): Different therapy schedules for Patient 2.

DownLoad: Full-Size Img PowerPoint

In Figure 6(A), we decreased the dosage given to Patient 1 to $20$, $19$, and $10$ percent of the original dosage, and depict the respective evolution with the green, red and cyan curves, respectively. At $19$, and $20$ percent of the original dosage, the proliferative cell mass steadily decreases during the therapy. While $20$ percent is enough to defeat the tumor before the therapy stops, for $19$ percent, a tiny amount of proliferative cells remain. These remaining cells then lead to a relapse of tumor growth, such that the tumor exceeds its largest volume even before our simulation time frame is over. At $10$ percent of the original dosage, we arrive at an equilibrium of the proliferative cells, in which the immunotherapy manages to control its growth but fails to eradicate it. Thus, as soon as the therapy is terminated, the proliferative tumor starts to grow again at the original rate. Finally, we note that even considering safety margins, $50$ percent of the dosage would have been more than enough to cure the cancer.

For Patient 2, an interesting question is how the changes in the therapy schedule might have affected the tumor growth. In Figure 6(B), in addition to the original therapy (blue curve), we have added a more aggressive therapy (red curve) where the patient is given the same dosage continuously every three weeks. In addition, we added a scenario where the switch to the six-weekly cycles was performed, but no therapy break occurred (green curve). We note that extending the therapy at the given point had close to no effect on the tumor growth. Even the more aggressive therapy only managed to slow down the growth dynamics by roughly four months before the tumor relapsed. It appears that, according to our model, the given tumor cannot be cured with the used drug, but the therapy only manages to prolong the patient's life by a certain time. Even significant increases in the dosage are not enough to eliminate the cancer, since the Hill–Langmuir equation which models the impact of the immunotherapy, approaches a stationary value for $\phi_\tau \to \infty$. With a continuous three-weekly drug administration scheme, further simulations show that increasing the dosage by a factor of 2 increases the time until the explosive growth occurs by 8 months. A factor of 10 could give the patient an additional year, before the tumor growth continues.

4.3.   Discussion of the patient-specific parameters

We conclude with a short discussion of the parameter values. Note that the dosage $d_\tau$, the drug's molar mass $M_\tau$ and the serum's half-life time $t_{{}^{1}\!\!\diagup\!\!{}_{2}\;}$ are known a priori from the given therapy, as well as the patient's weight $\omega$.

The parameters $\lambda^{ \rm{pro}}_P$, $\lambda$, $\lambda_{PN}$ and $\lambda^{ \rm{eff}}_\tau$ were assumed to be patient-specific and vary in our examples over a wide range. Concerning the choice of varying patient-specific parameters, we tried to limit ourselves directly connected to the physiological properties of the tumor. These are the tumor growth ($\lambda$, $\lambda^{\text{pro}}_P$), the efficiency of the immunotherapy ($\lambda^{\text{eff}}_\tau$) on the tumor cells, and the cell's reaction to malnourishment ($\lambda_{PN}$). In particular, $\lambda$ has a strong influence if the tumor stays localized, $\lambda = 2$, or spreads over the entire domain, $\lambda = 1$, leading to explosive growth later on. The parameters $\lambda^{ \rm{pro}}_P$ and $\lambda^{ \rm{eff}}_\tau$ have to counterbalance each other. Especially for Patient 2, where both growth and therapy are always active, they are difficult to determine. We remark, that the parameter $\lambda$ has a huge influence on the order of magnitude of $\lambda^{ \rm{pro}}_P$, and indirectly also on $\lambda^{ \rm{eff}}_\tau$. If we disable the immunotherapy, Patient 2 attains a volume of 1 cm³ at day 620 and doubles that volume at day 726 after 106 days, which is in a typical range for NSCLC. Therefore, even though growth parameters for both patients are very different, the growth itself is similar. We tried to keep those parameters related to the interfacial width of the Cahn-Hilliard equation fixed. Furthermore, we assume that the parameters of our nutrient model can somehow be calibrated patient-independently given enough data. Nevertheless, we are confident that better coefficients than the one presented can be found with more data. Especially the choice of $\eta_{iv}$ is primarily due to computational convenience.

To further assess the individual volatility of our parameters on the tumor, we conduct a local sensitivity analysis similar to the approach in [25,Sec. 4.1]. We define the local sensitivity indices of a scalar field $f$ dependent on the parameter $\alpha$ as

where $\bar{Q}(f)$ denotes the time averaged mass of the field $f$ and $g_{f, \alpha}$ its relative change up to first order when varying $\alpha$ by 1 percent. As scalar fields $f$ we consider the tumor $\phi_T$, the proliferative cells $\phi_P$, and the necrotic cell species $\phi_N$, as well as the visible tumor $V_{vis}$. For both patients, we set $T$ to the end of the respective immunotherapies. Numerically, we approximate $\partial_\alpha \bar{Q}(f)$ with a central finite difference with step size $h = \alpha \cdot 10^{-3}$ and evaluate the integral over time of $\bar{Q}$ with the trapezoidal rule.

The results for our two patients are depicted in Table 2: As expected, increasing the growth parameter $\lambda^{\text{pro}}_P$ increases the average cells mass while improving the therapy with $\lambda^{\text{eff}}_{\tau}$, or penalizing the growth of small cell volume fractions with $\lambda$ decreases the average mass. For these three calibrated parameters, $\lambda$ is the most sensitive parameter, followed by $\lambda^{\text{pro}}_P$, and finally $\lambda^{\text{eff}}_{\tau}$. It is noticeable that the cell mobility $c_m$ and growth factor $\lambda^{\text{pro}}_P$ have a much larger sensitivity for Patient 2 than for Patient 1. Here, both parameters increase how fast the tumor can spread over the entire lung domain, an effect which only appears for small $\lambda$ and is hence not relevant for Patient 1. In addition, we have more mixed populations of proliferative and healthy cells with Patient 2, which results in a larger impact of the nutrient consumption $\alpha_H$. As expected, a larger nutrient uptake of healthy cells results in smaller proliferative cell concentrations and increases the necrotic mass. On the other hand, $\lambda_{\text{PN}}$ and $\sigma_{\text{PN}}$ are more vital for Patient 1, since for Patient 2, the necrotic core mostly forms at the beginning due to the ongoing therapy. In contrast, for Patient 1, it develops over a longer period until the treatment takes effect.

Table 2.  Local parameter sensitivity analysis of Eq (4.1) for both patients and a subset of parameters. The value in brackets gives the absolute position compared to the other sensitivity values.

$\alpha$	$g_{\phi_T, \alpha}$		$g_{\phi_P, \alpha}$		$g_{\phi_N, \alpha}$		$g_{V_{vis}, \alpha}$
	Patient 1:
$\lambda^{\text{pro}}_P$	1.10E-02	(3)	9.02E-03	(2)	1.42E-02	(4)	1.02E-02	(5)
$\lambda^{\text{eff}}_\tau$	-3.55E-03	(7)	-5.08E-03	(6)	-9.57E-04	(10)	-5.16E-03	(7)
$\lambda$	-1.64E-02	(1)	-1.37E-02	(1)	-2.10E-02	(3)	-1.50E-02	(4)
$\lambda_{\text{PN}}$	1.77E-03	(10)	-2.20E-04	(10)	5.13E-03	(8)	2.32E-03	(10)
$\sigma_{\text{PN}}$	1.05E-02	(4)	-1.30E-03	(9)	3.03E-02	(2)	3.01E-02	(2)
$c_m$	-2.10E-04	(11)	-2.08E-04	(11)	-2.14E-04	(11)	-3.41E-04	(11)
$\varepsilon$	8.81E-03	(5)	6.78E-03	(3)	1.22E-02	(5)	7.66E-03	(6)
$c_\Psi$	-4.62E-03	(6)	-3.60E-03	(7)	-6.33E-03	(7)	-4.25E-03	(9)
$\kappa_i$	-2.87E-03	(8)	1.52E-03	(8)	-1.03E-02	(6)	-2.15E-02	(3)
$\alpha_H$	-2.51E-03	(9)	-6.53E-03	(4)	4.29E-03	(9)	-4.37E-03	(8)
$\alpha_P$	1.28E-02	(2)	-6.25E-03	(5)	4.50E-02	(1)	3.09E-02	(1)
	Patient 2:
$\lambda^{\text{pro}}_P$	4.13E-02	(2)	4.17E-02	(2)	6.19E-04	(11)	1.75E-03	(11)
$\lambda^{\text{eff}}_\tau$	-2.08E-02	(4)	-2.09E-02	(4)	-3.41E-03	(7)	-6.74E-03	(6)
$\lambda$	-1.82E-01	(1)	-1.84E-01	(1)	-1.71E-03	(8)	-5.02E-03	(7)
$\lambda_{\text{PN}}$	3.61E-04	(11)	3.28E-04	(11)	3.94E-03	(5)	2.17E-03	(10)
$\sigma_{\text{PN}}$	2.17E-03	(8)	1.97E-03	(8)	2.36E-02	(3)	1.33E-02	(3)
$c_m$	4.45E-03	(6)	4.50E-03	(6)	-8.46E-04	(10)	-2.79E-03	(8)
$\varepsilon$	1.52E-02	(5)	1.54E-02	(5)	-3.52E-03	(6)	-1.09E-02	(4)
$c_\Psi$	-3.15E-03	(7)	-3.18E-03	(7)	9.14E-04	(9)	2.54E-03	(9)
$\kappa_i$	-1.66E-03	(10)	-1.52E-03	(10)	-1.76E-02	(4)	-9.74E-03	(5)
$\alpha_H$	-3.64E-02	(3)	-3.71E-02	(3)	3.44E-02	(2)	1.71E-02	(2)
$\alpha_P$	2.14E-03	(9)	1.76E-03	(9)	4.35E-02	(1)	2.48E-02	(1)

---

 | Show Table

DownLoad: CSV

Since the nutrients control both the growth and decay of the proliferative cell populations, its diffusion factor $\kappa_i$ still has an impact on the tumor evolution. Similarly, since growth only happens at the interfaces, its width $\varepsilon$ and the scaling factor before the double well potential $c_\Psi$ also have a large influence.

4.4.   3D results

Finally, Figure 7 shows the tumor evolution for the full 3D model applied to Patient 1. On the top, a depiction of the lung mask in light gray, of the vasculature $\Gamma$ in red, a vertical and horizontal slice for $\phi_{\sigma, i}$, and the tumor in the lower right are shown. The proliferative tumor cells are colored blue, while the necrotic cells are depicted in yellow. A threshold of $0.3$ is chosen for both the contours of the proliferative and necrotic cells. For the given snapshot in time, the shell of proliferative cells is seen to be oriented in the direction of the nutrient gradient. On the right, the nutrient concentration at different times is shown. Already at day 10, a nutrient shortage close to the tumor is observed. This area grows with the tumor until it reaches its maximal size after 140 days, shortly before the immunotherapy starts. After 420 days, most of the proliferative cells consuming nutrients have disappeared, and the effects of the tumor on the nutrient concentration are no longer visible. The tumor evolution for certain points in time is shown. After 60 days, a small necrotic core has formed. It is not at the center of the tumor, but slightly shifted towards the nutrient-poor domain. Again, the proliferative cells further away from the nutrients decay until only a tiny proliferative shell remains, which points towards the nutrient source. Finally, only a few proliferative cells remain inside the tumor, mainly consisting of necrotic cells. At this point, the visible tumor volume no longer changes and the proliferative core decays further.

Figure 7.  3D tumor simulation: Top: Nutrient concentration $\phi_{\sigma, i}$ after 200 days have elapsed. Bottom: Tumor shapes consisting of proliferative (blue) and necrotic (orange) cells at different points in time.

DownLoad: Full-Size Img PowerPoint

We want to stress that even though both 1D and 3D produce quantitatively the same results in terms of tumor volumes, we expect that the 3D simulation will turn out superior if more knowledge about the tumor and the surrounding tissue is available and used for a refined nutrient model. Already for the given simplistic model, we see that the simulated 3D tumor is generally not spherically symmetric which violates the 1D assumptions leading to very crude approximations. This is even more the case for the metastases of Patient 2, where it is quite relevant for the treatment where the tumor reappears.

5.   Discussion

Critical microenvironments influencing tumor formation are generally not observable ^[78]. Deterministic modeling approaches, on the other hand, can aid in understanding the driving dynamics by simulating growth and decline for observable tumor sizes. It has been demonstrated here, that our phase-field model can qualitatively describe the tumor volume evolution in the observable window over time for NSCLC patients and can address different outcomes of immunotherapeutic treatment approaches. The simulation results of two clinical cases were explained in a biologically meaningful manner throughout the observation period with results that agree with multiple volumetric measurements at time points acquired during clinical routine examinations. This is considered as a preliminary proof of concept for the presented model. In future studies, where the model is calibrated over larger data sets, the parameter estimation should be addressed in particular. The parameter estimation is then expected to be more stable if prior parameter distributions for general patient cohorts can be identified.

We have shown that our model easily generalizes to a full simulation in 3D. Backdrawn data on RECIST measurements, which are clinically routine characteristics, are easily processable, but the presented model may contain significantly more clinically relevant information by providing the spatial structure of tumors under therapy influence. However, the 1D representation is potentially already very useful for a swift evaluation of qualitative treatment outcomes of different immunotherapeutic dosages.

The model simulation strongly depends on patient-specific parameters unknown a priori and inferred over the whole simulation time, including the clinical outcomes. However, if these parameters can be estimated quickly during treatment or are based on patient-specific clinical covariates, then this approach has prognostic potential. The long-term therapeutic success and lower but still successful dosages could possibly be determined. A preview of this potential is presented in our results, as the qualitative therapeutic outcome can be determined with alternative drug dosage regimens for the two patients presented.

This potential should be addressed in future work and be based on larger patient data sets. Successful early estimation of these parameters allows for an improved treatment schedule and a minimal drug dosage regimen to reduce adverse effects. In addition, it could also reduce data acquisition efforts by reducing the number of CT scans during patient treatment to a minimal amount needed for model control, re-calibration, and clinical check-ups. This has the potential to increase the patient's quality of life on an individual basis.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgements

The work of M. Fritz, A. Wagner and B. Wohlmuth was partially funded by the Deutsche Forschungsgemeinschaft (WO 671/11-1). P. Schlicke was partially funded through IGSSE/TUM Graduate School. M. Fritz is partially supported by the State of Upper Austria. The support of J. Tinsley Oden by the U.S. Dept. of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program, under Award DE-960009286 is gratefully acknowledged.

Conflict of interest

All authors declare no conflicts of interest in this paper.

Code and data availability

The data and code used for the simulations are available at https://github.com/wagnandr/immunotherapy-lung-cancer.