PaperlayerView exact revision on arXiv

Source-linked AI summary

Cahn-Hilliard phase-field modeling of tumor growth via locally adaptive isogeometric analysis with THB-splines

Dhiraj S. Bombarde, Carlotta Giannelli, Alessandro Reali, Guillermo Lorenzo

arXiv:2609.08479v2math.NAcs.NA

TL;DR

Predicting tumor dynamics under physiologically relevant conditions remains challenging, and the higher-order Cahn–Hilliard formulation imposes demanding continuity requirements on primal discretizations. This paper combines the model with locally adaptive isogeometric analysis using THB-splines. Its benchmarks reproduce established tumor morphologies, and the approach is applied to an MRI-based patient-specific breast geometry.

  • Problem

    Predicting tumor dynamics under physiologically relevant conditions remains challenging, while the fourth-order Cahn–Hilliard operator requires at least C1 continuity for a well-defined primal variational formulation.

    Source: §1 · pp. 1–2

  • Method

    The paper solves a nutrient-driven Cahn–Hilliard tumor growth model in primal form using isogeometric analysis with locally adaptive THB-splines.

    Source: §5 · p. 23

  • Results

    The approach reproduces spheroidal and finger-like tumor morphologies in 2D and 3D benchmarks and is applied to an MRI-reconstructed patient-specific breast geometry.

    Source: §1 · p. 3, §4.3 · p. 23

  • Takeaways & Limitations

    The results demonstrate the potential of the Cahn–Hilliard phase-field model and locally adaptive IGA framework for tumor-growth simulations in clinically relevant breast geometries.

    Source: §1 · p. 3, §4.3 · p. 23

    Without calibration against longitudinal MRI data, the framework’s predictions remain qualitative rather than quantitative and patient-specific.

    Source: §5 · p. 23

Abstract

from arXiv · show

Predicting tumor dynamics in biological systems under physiologically relevant conditions using mathematical and computational models remains a challenging problem. Continuum models based on phase-field (diffuse-interface) formulations have proven to be an effective modeling strategy to govern tumor dynamics and interactions of multiple species. Within this framework, the present work investigates a tumor growth model based on the Cahn-Hilliard (CH) equation. The formulation involves a fourth-order differential operator that imposes higher continuity requirement on approximation spaces for a well-defined primal variational formulation. To address this challenge, we use isogeometric analysis (IGA), which inherently satisfies this requirement through spline-based basis functions and eliminates the need for mixed or auxiliary-variable approaches commonly used in standard finite element discretizations. Additionally, a locally adaptive IGA scheme with truncated hierarchical B-splines (THB-splines) is used to reduce computational cost while maintaining accuracy. The model is first evaluated on standard benchmark cases and then applied to an organ-scale, patient-specific geometric model of the breast reconstructed from magnetic resonance imaging (MRI) data. Our results show that the model reproduces known tumor morphologies, ranging from a spheroidal pattern to fingered growth. A series of numerical experiments further shows the diversity of tumor dynamics produced by different model parameter choices. Our findings demonstrate the predictive potential of the CH-based phase-field tumor growth model integrated with a locally adaptive IGA framework.

1. Introduction

The paper targets the difficulty of predicting tumor dynamics under physiologically relevant conditions and presents a primal Cahn–Hilliard model solved with locally adaptive isogeometric analysis. It evaluates the approach on 2D and 3D benchmarks and a patient-specific breast geometry.

1. Introduction

  • Predicting tumor dynamics under physiologically relevant conditions remains challenging, while personalized disease management increases the need to anticipate long-term disease evolution.Continuum phase-field formulations represent interacting species and evolving tissue interfaces within a unified description.Source: §1 · pp. 1–2
  • Cahn–Hilliard models maintain physically realistic diffuse interfaces over time, prevent unphysical mass changes, and can couple consistently with nutrient diffusion and other processes.Their fourth-order operators require at least C1 continuity for a well-defined primal variational formulation, a constraint that complicates finite-element discretization, especially in 3D.Source: §1 · pp. 2–3
  • The proposed approach solves the Cahn–Hilliard tumor model in primal form with isogeometric analysis, whose spline basis functions satisfy the higher-continuity requirement.This avoids the mixed or auxiliary-variable formulation commonly used for finite-element solutions of higher-order operators.Source: §1 · pp. 2–3
  • Locally adaptive THB-splines refine steep phase-field interfaces and coarsen less active regions to balance accuracy with degrees of freedom and computational effort.Their reduced support also improves system-matrix sparsity compared with standard hierarchical B-splines.Source: §1 · p. 3
  • The model reproduces spheroidal and finger-like tumor morphologies in 2D and 3D benchmarks and is applied to an MRI-reconstructed patient-specific breast geometry.The authors describe this as the first higher-order Cahn–Hilliard model of breast cancer growth in a realistic, organ-scale 3D setting.Source: §1 · p. 3

2. Model formulation

This section formulates a nutrient-driven Cahn–Hilliard phase-field model for tumor growth and gives the weak forms used for its numerical discretization. The model tracks tumor and nutrient fields, with stated assumptions on nutrient transport, boundary conditions, and omitted mixture velocity.

2. Model formulation

  • The model couples tumor growth through proliferation, apoptosis, and chemotaxis with nutrient dynamics, and presents governing equations, boundary and initial conditions, and their weak formulation.The weak formulation provides the basis for the numerical discretization in later sections.Source: §2 · p. 3

2.1. CH-based tumor growth model

  • The model uses a tumor phase field ϕ governed by a Cahn–Hilliard-type equation and a nutrient concentration σ governed by a reaction–diffusion equation.The two fields describe a mixture of healthy and tumor tissue.Source: §2.1 · p. 3
  • The order parameter ϕ is the difference between the tumor and healthy tissue volume fractions, with ϕ = −1 and ϕ = 1 denoting healthy and tumor tissue, respectively.The mass-balance expression assumes constant densities and no excess volume due to mixing.Source: §2.1 · p. 4
  • The phase-field equation uses a chemical potential derived from free energy, with nutrient-induced energy effects and chemotaxis represented in the energy terms.The diffuse-interface thickness is controlled by ϵ, and the source term accounts for nutrient-dependent proliferation and apoptosis.Source: §2.1 · p. 4
  • The source term follows linear kinetics, with nonnegative constants P and A associated with proliferation and apoptosis, respectively.The nutrient sink is assumed linear in concentration and phase, and nutrient dynamics are quasi-static because diffusion is much faster than cell proliferation.Source: §2.1 · p. 4
  • Nutrient transport is represented by a diffusive flux and a consumption sink; the model assumes constant nutrient diffusivity and prescribes boundary nutrient values to represent supply from well-perfused healthy tissue.Homogeneous Neumann conditions are used for ϕ and its chemical potential µ.Source: §2.1 · pp. 4–5
  • The volume-averaged mixture velocity considered in the referenced model is excluded to keep the present formulation simple.The paper identifies its inclusion as a natural extension.Source: §2.1 · p. 5

2.2. Variational formulation of the coupled system

  • The coupled system is written in a Bubnov–Galerkin weak form, with a staggered phase-field solve for a given nutrient field and weak enforcement of the phase-field boundary condition using Nitsche’s method.The nutrient equation is also formulated variationally; spatial discretization uses Galerkin IGA and time integration uses the generalized-α scheme.Source: §2.2 · p. 5

3. Adaptive IGA with THB-splines

The section develops a locally refinable THB-spline IGA framework for the fourth-order Cahn–Hilliard model, covering hierarchical meshes, error-guided adaptation, spatial discretization, and time integration.

3. Adaptive IGA with THB-splines

  • The adaptive IGA framework recalls spline basics and combines locally refined THB-spline meshes with refinement and coarsening strategies that preserve suitable grading.Its hierarchical mesh comprises active elements drawn from nested refinement levels.Source: §3 · p. 5, §3.2 · p. 6

3.1. B-splines

  • Multivariate B-spline bases are formed as tensor products of univariate B-splines defined over knot vectors.The univariate basis is defined recursively by the Cox–de Boor formula.Source: §3.1 · p. 6

3.2. THB-splines

  • THB-splines truncate hierarchical basis functions to reduce their support, improving matrix sparsity and enabling more localized refinement with fewer degrees of freedom than HB-splines.The hierarchical construction uses nested domains and spline spaces across refinement levels.Source: §3.2 · pp. 6–7

3.3. Error estimator

  • A gradient-based element error estimator and tolerances α and β guide element marking for refinement or coarsening using the greatest appearing error utilization strategy.The local error uses the squared L2 norm of the phase-field gradient over an element and its volume.Source: §3.3 · p. 7

3.4. Admissible THB-spline refinement and coarsening

  • An admissible mesh of class m limits the basis functions nonzero on any element to at most m consecutive levels, keeping their overlap independent of the total refinement levels.Coarsening is conservative: all children must be marked, and the coarsening neighborhood must be empty to reactivate an element.Source: §3.4 · pp. 7–8

3.5. Spatial discretization

  • THB-spline spaces with at least C1 continuity directly discretize the fourth-order Cahn–Hilliard operator, while the nutrient field uses the same space for consistency and stability.This avoids introducing auxiliary variables for the phase-field equation.Source: §3.5 · p. 9

3.6. Time-stepping scheme

  • The generalized-α method advances the nonlinear phase-field system using Newton–Raphson iterations, after which the nutrient concentration is computed from a linear system.The nutrient update follows convergence of the phase-field solution at each time step.Source: §3.6 · p. 10

3.7. Adaptive refinement and coarsening

  • Adaptive meshes are initialized from an error estimator and periodically refined and coarsened; the updated solution is transferred between spaces by L2 projection.Refinement captures evolving features, while coarsening removes superfluous resolution to reduce computational overhead.Source: §3.7 · pp. 10–11

4. Numerical studies

The numerical studies test adaptive THB-spline IGA on 2D and 3D benchmarks, then apply it to an MRI-derived breast geometry. They reproduce multiple tumor morphologies, while showing how parameter choices, initial conditions, and the patient-specific setting shape the results.

4.1. Tumor growth in a square domain

  • Adaptive meshes with ℓ=6 and ℓ=5 yield nearly identical L2 and H1 errors; ℓ=6 uses fewer DOFs and similar or marginally lower runtime overall.Source: Figure 7 · p. 14, Figure 8 · p. 14, §4.1 · pp. 15–16
  • The selected square-domain mesh reproduces finger-like protrusions in close agreement with the overkill reference and results reported in.A p=3 mesh with 1024 elements per direction is adequate for the deepest interface refinement.Source: §4.1 · pp. 15–16, Figure 10 · p. 16
  • In the square-domain tests, increasing proliferation from P=0.1 to 5 produces smooth spheroidal growth, while other parameter adjustments can produce elongated shapes.Source: §4.1 · pp. 16–17, Figure 11 · p. 17
  • With parameters fixed, minimal perturbations to the square-domain initial profile can produce entirely different growth morphologies.Source: §4.1 · p. 17, Figure 12 · p. 18

4.2. Tumor growth in cubic domain

  • In the cubic domain, an ellipsoidal seed develops an umbrellalike morphology, combining fingered vertical-plane and circular horizontal-plane growth, with qualitative validation against.Source: §4.2 · p. 18, Figure 13 · p. 19
  • In the cubic case, increasing proliferation from 0.1 to 5 shifts growth toward a spheroidal, volumetric morphology and steepens nutrient depletion in the tumor interior.The authors attribute the steeper depletion profile to increased consumption at higher proliferation.Source: §4.2 · p. 18, Figure 14 · p. 20

4.3. Tumor growth in patient-specific organ-scale breast model

  • The MRI-derived breast model assembles 18 traced slices into an IGA-compatible volume that closely agrees with a 3D Slicer reconstruction; the authors describe this as their first higher-order CH study in a realistic organ-scale 3D setting.Source: §4.3 · pp. 19, 21, Figure 16 · p. 21
  • In the patient-specific breast, an elongated boundary-adjacent tumor undergoes migration-dominated evolution with limited morphological change.The authors note that more realistic dynamics could be achieved by calibrating model parameters to longitudinal MRI data.Source: Figure 18 · p. 22, §4.3 · p. 23

5. Conclusion

The work applies primal Cahn–Hilliard tumor modeling with adaptive THB-spline IGA to benchmark morphologies and a patient-specific breast geometry. The authors identify longitudinal-MRI calibration as a future step toward quantitative patient-specific predictions.

5. Conclusion

  • Adaptive THB-spline IGA refines the evolving diffuse interface and coarsens low-activity regions to solve the nutrient-driven Cahn–Hilliard model in primal form.Numerical experiments establish mesh parameters, and 2D and 3D benchmarks examine complex morphologies and sensitivity to model parameters and initial conditions.Source: §5 · p. 23
  • The framework reproduces a broad range of clinically relevant tumor morphologies in benchmarks and is applicable to organ-scale simulations on patient-specific breast geometry reconstructed from MRI.The authors describe this as demonstrating predictive potential in patient-specific settings, not as quantitative clinical validation.Source: §5 · p. 23
  • Longitudinal MRI calibration is identified as a future requirement for moving from qualitative morphology prediction to quantitative, patient-specific growth modeling.The authors also propose mechanical feedback from surrounding breast tissue and treatment-response modeling as possible extensions.Source: §5 · p. 23

Appendix A. Implementation workflow and pseudocode

Algorithm A1 outlines the computational workflow for the phase-field–nutrient system, combining nonlinear solution and time integration with adaptive mesh refinement and coarsening.

Appendix A. Implementation workflow and pseudocode

  • Algorithm A1 traces the workflow from initialization through nonlinear solution and time integration to adaptive mesh refinement and coarsening.The listed scheme uses the NR method for nonlinear solution and the generalised-α scheme for time integration.Source: Appendix A · p. 24
  • The workflow initializes the phase field and nutrient state, then iterates the phase-field solve and updates the nutrient solution at each time step.The pseudocode includes a nonlinear iteration loop and a nutrient-equation solve.Source: Appendix A · p. 24
  • At adaptation steps, the algorithm evaluates an element indicator, marks elements for refinement or coarsening, enforces admissibility class m, and transfers the solution.The coarsening step uses an L2 projection to obtain an updated phase-field solution.Source: Appendix A · p. 24

Appendix B. Extended numerical studies

Appendix B varies model parameters and initial conditions, tests quadratic THB-splines, and considers an interior tumor in patient-specific breast geometry. The cases yield distinct morphologies, while quadratic bases do not attain the accuracy of cubic or quartic counterparts.

Appendix B. Extended numerical studies

  • With E = 5, λ = 0.002, and M = 5, the same initial condition produces an elliptical tumor; increasing λ broadens the diffuse interface.The adopted mesh resolution is described as marginally conservative for this configuration.Source: Appendix B · p. 25, Figure B.1 · p. 25
  • A cosine perturbation with six-fold angular modulation of amplitude 1/40 breaks initial radial symmetry and leads to a distinct fingering pattern.The interface instabilities are amplified by the underlying dynamics.Source: Appendix B · p. 25, Figure B.2 · p. 26
  • Quadratic B-splines cannot match the accuracy of cubic or quartic counterparts for the tested interfacial features, even at the best attainable quadratic solution quality.The authors frame quadratic bases as an option when a moderate reduction in accuracy is acceptable.Source: Appendix B · pp. 25–27
  • In the breast geometry, an interior ellipsoidal tumor experiences less influence from nutrient-rich boundaries, making growth-driven changes more apparent and migration no longer dominant.The initial ellipsoid has two semi-axes of 0.11 and a third of 0.13.Source: Appendix B · p. 27, Figure B.4 · p. 27
Loading 2609.08479v2…