Computational Modeling of Patient-Specific Healing and Deformation Outcomes Following Breast-Conserving Surgery Based on MRI Data.

Introduction
Breast cancer is among the most common cancers in women, with roughly 1 in 8 women diagnosed during their lifetime and over 367,000 new cases reported annually in the USA [1]. Increased awareness, routine screenings that enable early detection, and advancements in treatment have significantly improved survival rates, with localized breast cancer now having a 5-year survival rate exceeding 99% [1]. With these improved outcomes, clinical focus is shifting toward long-term considerations, including recurrence prevention and enhancement of quality of life.
Surgical treatments are associated with the lowest recurrence rates [2, 3] and are typically performed via breast-conserving surgery (BCS; otherwise known as lumpectomy) or mastectomy (complete removal of the breast). Choosing between these options is often complex and emotionally challenging, involving personal preferences, patient and cancer characteristics, and anticipated post-surgical outcomes. In recent years, BCS has emerged as the preferred standard of care for early-stage breast cancer, offering recurrence and survival rates comparable to mastectomy while carrying a lower risk of complications [4–8]. The primary goal of BCS is to completely remove the cancer while preserving as much healthy breast tissue and appearance as possible. However, the resulting surgical cavity initiates a complex healing process, resulting in variable levels of tissue contraction, scar tissue formation, and breast deformation (i.e., cosmetic defects, including dents, distortions, and asymmetries between breasts). Achieving a favorable cosmetic outcome is a critical component of surgical decision-making, as it is strongly associated with psychologic recovery and long-term quality of life [9, 10]. Despite its importance, predicting oncologic, healing, and cosmetic outcomes remains challenging due to the complex nature of breast tissue repair and inter-patient variability. Factors such as breast and tumor geometry, tumor location, and breast density are known to influence healing, yet their contributions remain poorly understood. A deeper understanding of how these patient-specific characteristics affect BCS outcomes is needed to support more personalized treatment plans, reduce post-surgical complications and secondary procedures (i.e., re-excision or reconstruction), and improve overall satisfaction and quality of life [4].
Fundamentally, the healing process of the post-BCS surgical cavity follows the classic overlapping phases of reparative wound healing: hemostasis, inflammation, proliferation (or granulation), and remodeling, as demonstrated in a porcine lumpectomy model [11]. Healing begins immediately after the cavity is formed, with the surgical dead space typically filling with blood (hematoma) and/or serous fluid (seroma) [11]. A provisional fibrin matrix then forms, which lacks mechanical integrity, allowing the cavity and surrounding tissue to contract while promoting inflammation and a cytokine gradient [11]. This gradient promotes angiogenesis, fibroblast proliferation and infiltration, and collagen deposition [11]. As fibroblasts differentiate into myofibroblasts, they further drive tissue contraction and fibrotic remodeling, resulting in stiff scar tissue that can lead to pain, breast deformities, and altered breast consistency, all of which can negatively impact a patient’s psychologic well-being [12].
Over the past three decades, mechanistic computational models have been developed to depict the wound healing process, primarily focusing on cutaneous (skin) repair. These models have evolved to better capture some of the multi-scale and multi-phase complexities of the healing process [13, 14]. More recently, computational models have been applied to BCS, aiming to predict post-surgical healing and cosmetic outcomes [15–17]. However, these models lack comprehensive calibration against experimental and clinical breast healing data. Furthermore, current models fail to capture the complex interactions between cellular mechanobiological activity, extracellular matrix (ECM) deposition and remodeling, and tissue contraction, all which are critical factors for accurately modeling the long-term contracture of breast tissue and overall breast deformation.
As an initial step to address this gap, we previously developed and calibrated a computational mechanobiological model of breast cavity healing following BCS using porcine lumpectomy data and published clinical data [18]. As detailed in Sect. 2.2, the finite element model integrates biologic, microstructural, and mechanical dynamics, including fibroblast infiltration, collagen remodeling, and ECM stiffening, capturing key drivers of breast tissue contraction and deformation. The model also couples nonlinear tissue mechanics with biologic fields (i.e., fibroblast density and inflammation signals) described through reaction–diffusion equations. Pro-inflammatory cytokines drive fibroblast infiltration and proliferation, which in turn promote local collagen deposition, ultimately influencing tissue stiffness and active stresses contributing to contraction. Our previous simulations used a generalized breast geometry [18]. Thus, in our original work, while we were able to calibrate the dynamics of the wound healing model, we were unable to investigate the effect of patient-specific factors, namely breast shape, density, and tumor location and size.
In parallel with computational efforts, several clinical studies have evaluated how patient-specific characteristics influence BCS outcomes using statistical analyses of patient-reported data. For example, surgical decision-making models have correlated breast and tumor characteristics with post-surgical esthetics and quality of life outcomes [19–21]. These studies show that factors such as tumor-to-breast volume ratio and tumor location can significantly affect cosmetic results and patient well-being, offering potential guidance on whether a patient may benefit more from BCS or mastectomy. Additional work using questionnaires and other surveys has revealed associations between variables such as lower breast density, higher cavity-to-breast volume percentage (CBVP), increased age, elevated body mass index (BMI), breast irradiation, and adjuvant therapy with poorer cosmetic outcomes and reduced patient satisfaction [22–27]. However, while these studies identify key trends, their reliance on questionnaire-based feedback limits insight into the underlying mechanobiological processes driving individualized healing outcomes.
In this study, we build upon our previously developed and experimentally calibrated model of post-BCS breast cavity healing [18], integrating a large dataset of preoperative magnetic resonance imaging (MRI) scans from the Cancer Imaging Archive [28] to define how individual factors, such as breast geometry, breast density, cavity volume, and cavity location, impact breast healing and deformation outcomes. As detailed in Sect. 2.3, simulation results were further used to train machine learning surrogate models, enabling rapid evaluation of the effects of patient-specific features on healing dynamics and cosmetic outcomes. These patient-specific simulations and machine learning (ML) metamodels have the potential to support breast surgeons and patients in making more informed surgical decisions and developing personalized treatment plans, with the goal of enhancing outcome predictability, minimizing complications, and improving patient satisfaction and quality of life.

Materials and Methods

Preoperative MRI Processing
Patient-specific data were obtained from preoperative MRI scans of breast cancer patients, sourced from The Cancer Imaging Archive, with data collected by Duke University Hospital [28]. A total of 103 patients were randomly selected from the database, excluding those with bilateral breast cancer or multiple tumors in the same breast. As displayed in Fig. 1, each MRI was processed using 3D Slicer, with tissue segmentation performed by applying intensity thresholding to estimate the volumes of adipose tissue, fibroglandular tissue, and tumor tissue [29]. Breast tissue regions were defined relative to the chest wall using MRI intensity and anatomic landmarks to distinguish boundaries from surrounding thoracic structures. Adipose and fibroglandular tissue were then segmented in 3D Slicer by applying MRI intensity thresholding, with their respective volumes quantified. Breast density for each patient was calculated as the percentage of fibroglandular tissue relative to the total breast volume. Three-dimensional breast geometries were also generated from the MRI data and meshed in 3D Slicer. Tumors were modeled as ellipsoids using dimensions provided in the repository’s datasheets. A summary of the patient and tumor characteristics is provided in Table 1.
Since MRIs are acquired with patients lying in a prone position, the breasts are deformed under gravity, which does not reflect their natural upright shape and appearance. To correct for this, finite element simulations were performed to model the deformation of the breast and tumor as they transition from the prone to the upright position (Fig. 1). Both the breast and tumor were modeled as nearly incompressible, neo-Hookean hyperelastic materials. Material properties of the breast tissue were derived from each patient’s breast density using the rule of mixtures. Young’s modulus values of 1.5 kPa for adipose and 6 kPa for fibroglandular tissue were chosen, consistent with the work of Eder et al. (2014), who tested a range of reported material properties in finite element simulations of prone-to-upright breast deformation against 3D surface scans of patients [31]. The post-surgical breast cavity was modeled using the tumor geometry and location. Material properties assigned to the cavity reflected serous fluid accumulation immediately after BCS, with a shear modulus of 0 kPa and a bulk modulus of 300 kPa to simulate its near-incompressibility.
The cavity volume for the 103 MRI-derived patient-specific geometries was based on standard surgical guidelines for early-stage invasive breast cancer, which recommend "no ink on tumor" and no additional negative margin width [32]. However, recognizing the variability in surgical techniques across different surgeons and patient cases, the extent of tissue removal can vary significantly. To account for this, an additional 60 geometries were generated using the same preoperative MRIs, with a randomized uniform negative cavity margin ranging from 0 mm and 5 mm, reflecting the range commonly used by surgeons [32, 33].

Computational Model Framework
The computational mechanobiological model is a custom finite element solver executed in C++, with the code repository link included at the end of the manuscript. The framework was initially developed to simulate cutaneous wound healing [34–36] and was subsequently adapted and calibrated to be specific to the breast [18]. A general overview of the model is provided below, with a more detailed description available in our previous publications [18, 35, 36]. Model parameters used in this study are consistent with those described in Harbin et al. (2023) [18]. In that work, 5 parameters governing fibroblast and collagen density were calibrated using experimental porcine lumpectomy data, while 4 parameters relating to active stress and plastic deformation were calibrated using human clinical cavity contraction data and subsequently evaluated through sensitivity analyses. Remaining parameters were informed from previously published experimental data and literature values [18]. The full list and description of model parameters, along with a brief explanation of how they were determined in our previous work [18], are provided in Supplementary Table S1.

Kinematics
The reference breast geometries are described with material coordinates . The time-dependent configuration, , is obtained through deformation mapping, , described as . Fibroblast density, cytokine concentration, and collagen density are defined as , respectively, with the collagen matrix described through fiber dispersion and preferred fiber orientation . Local geometries changes are captured through the deformation gradient, , which can be split to capture both elastic and plastic deformation . The plastic deformation tensor can be described through an orthonormal basis []Evolution of the various fields and parameter values are discussed in detail below and in the Supplemental information.

Biochemical Fields
Fibroblast density and cytokine concentration are represented as reaction–diffusion partial differential equationsThe fibroblast diffusion is dependent on both collagen density and cytokine concentrationwhile the cytokine diffusion coefficient is assumed constant. In this expression, represents fibroblast velocity as a function of collagen density and acts as a scaling constant that sets the magnitude of the collagen-dependent contribution. These parameters were calibrated in our previous work using spatiotemporal fibroblast density data from histologic analysis of the porcine lumpectomy study [11, 18]. The fibroblast and cytokine source terms are defined asCytokine-increased fibroblast proliferation, represented through in the fibroblast source term, was also calibrated to spatiotemporal fibroblast density data from the porcine lumpectomy study [11, 18]. The term , included in both source terms, represents a mechanosensing activation function:

Mechanical Equilibrium and Stress
In the absence of body forces, the balance of linear momentum is reduced towith appropriate boundary conditions. In the finite element simulations, the chest wall and contralateral breast were constrained with Dirichlet boundary conditions, while the exterior surface of the ipsilateral breast was modeled with a free boundary. The total stress comprises both passive and active contributionsdictated by the passive material response and active extracellular matrix (ECM) contraction and reorientation directed by fibroblasts and myofibroblasts. The material behavior is described through the hyperelastic strain energy functiontaking into account the isotropy of the non-collagenous matrix of the ECM and the possible anisotropy of the collagen fibers. It is parameterized by , with the neo-Hookean contributions, and , determined through each patient’s unique breast density using the rule of mixtures. Several studies were used to inform the material properties of adipose and fibroglandular tissue, as we estimate their Young’s modulus to be 10 kPa and 40 kPa, respectively [37–45].
The passive material response is also dependent on microstructure fields and elastic invariants of the deformation . These depict the elastic volume change, , the trace of the elastic right Cauchy Green tensor, , and the deformation in the preferred fiber direction, , with representing the deformed fiber orientation. The microstructure parameter , volume fraction of collagen, scales the strain energy, whereas and control anisotropy. Here, we consider initially which describes isotropic behavior regardless of the value of . As noted below, over time, anisotropy can develop as the fiber orientation becomes aligned with the direction of maximum deformation and the dispersion evolves to match the anisotropy of the deformation.
The passive stress is determined through the strain energy function byThe active stress features couplings between the fibroblast density, cytokine concentration, collagen density, and preferred fiber orientation:with , . In previous work, parameters and , representing the contractile force of fibroblasts and myofibroblasts, along with , the saturation of mechanical force by collagen, were optimized against longitudinal clinical data describing breast cavity contraction following lumpectomy [18].

Tissue Remodeling
Local remodeling of the ECM is encoded through a series of ordinary differential equations describing the microstructural fields . Collagen deposition is a function of both fibroblast density and cytokine concentrationCollagen production, , and cytokine-increased collagen production, , were previously calibrated to spatiotemporal collagen density data obtained from histologic analysis of the porcine lumpectomy study [11, 18]. The change in plastic deformation occurs independently in the established orthonormal basiswith . It is dependent on the positive term of the collagen rate change, , which signifies the new collagen deposition rate. The Macaulay brackets specify that plastic deformation only occurs beyond a set threshold deformation . Similar to the parameters relating to active stress, the rate of plastic deformation, , was also previously optimized to fit time-dependent human clinical data of breast cavity contraction [18].
The change in preferred collagen fiber orientationallows for the reorientation of the principal fiber direction to the direction of the maximum principal stretch, with being the largest eigenvalue and corresponding eigenvector. The change in fiber dispersionis dependent on the ratio of the first two eigenvalues, and . Collagen fibers were initialized at , representing isotropy, corresponding to an isotropic distribution with no preferred orientation. Anisotropy corresponds to lower dispersion values, indicating increasing alignment of fibers along a preferred direction.

Multi-Fidelity Gaussian Process Surrogates

Multi-Fidelity Gaussian Process Regression Overview
The finite element model is computationally expensive, making tasks such as sensitivity analyses impractical for large-scale patient-specific simulations. Furthermore, even though the random sampling of patients from the MRI database captures the breadth of possible patient-specific characteristics, for a new patient, a full finite element simulation would still be needed. To overcome these limitations, we employed multi-fidelity Gaussian process (GP) regression to develop a computationally efficient surrogate model.
The multi-fidelity GP regression approach is adapted from Perdikaris et al. (2017) [46], and an overview is shown in Fig. 2. The framework builds on standard GP regression with an autoregressive stochastic formulation and extends it to incorporate multiple datasets from multiple fidelity levels [47–49]. It is classified as a deep GP, in which the posterior distribution from a lower-fidelity level is used as a training input for the higher-fidelity level. This hierarchical structure allows nonlinear correlations between fidelities to be captured effectively. By leveraging these correlations, the model combines low- and high-fidelity simulation data to enhance prediction accuracy while significantly reducing computational cost.
The approach assumes that there are multiple datasets of varying fidelity expressed as , with t denoting increasing fidelity levels (t = 1, 2). represents the set of M total inputs (), with the corresponding outputs. To improve numerical stability, small Gaussian random noise with zero mean and variance is added to the output [50]. The lowest fidelity level can be modeled using standard GP regression, with the function assigned to a zero mean GP prior and k being a squared exponential covariance function assigned with hyperparameters . The hyperparameters define the covariance structure in terms of output variance, input length scales, and noise variance and were optimized by maximizing the marginal log-likelihood expressed aswhere is a symmetric, positive-definite matrix . The Gaussian process framework and hyperparameter optimization were carried out in Python using GPy, with convergence determined by tolerance on log-likelihood change or a maximum iteration limit [51]. Through Bayes’ rule, the posterior distribution can be found,and used for new predictions for output given the input . The predictive mean and variance areandrespectively, with .
For the subsequent fidelity level, an autoregressive formulation is utilized that adds an extra dimension to the training input for the next fidelity level in the form of the GP posterior from the previous fidelity levelwhere . Such an approach decouples the problem into standard GP regression schemes at each fidelity level, thus capturing nonlinear correlations between fidelity levels without the limitation of computational cost. Similarly, the covariance function is also a function of the posterior of the previous fidelity levelwith , , and being squared exponential covariance functions along with their associated hyperparameters , , and . The hyperparameters can again be learned through the maximization of the marginal log-likelihood, displayed in Eq. 17.
As stated above, the predictive posterior distribution for the lowest fidelity level is Gaussian with a corresponding mean and variance (Eqs. 19 and 20). The sequential fidelity level also yields a Gaussian posterior, as it is dependent on the Gaussian prediction from the first level, and is given byMonte Carlo integration is used to calculate the posterior predictive mean and variance for all fidelity levels following the lowest fidelity level [46].

Multi-Fidelity Gaussian Process Surrogate Training
In addition to the 163 high-fidelity simulations based on patient-specific breast geometries, 273 low-fidelity simulations were performed using a simplified generalized breast geometry with a coarser mesh. The low-fidelity meshes utilized approximately 70% fewer elements than the high-fidelity meshes, significantly reducing computational cost. Although less accurate than the high-fidelity simulations, the low-fidelity simulations captured the overall breast healing and deformation response, which was leveraged through correlations with the high-fidelity data to improve surrogate model accuracy. In the low-fidelity simulations, the breast and cavity geometries were approximated as a hemisphere and a sphere, respectively, as shown in the example in Fig. 3. Four patient-specific characteristics were used as input parameters for the surrogate models: breast density, cavity volume, breast volume, and cavity depth (defined as the shortest distance between the breast surface and the center of the cavity) (Fig. 3). These features were selected for their relevance to breast tissue mechanical properties, breast geometry, and cavity geometry and location. Input values were generated using Latin hypercube sampling [52] to ensure broad coverage of the characteristic space representative of the MRI-derived patient population (Table 1).
Two multi-fidelity GP surrogates were constructed as shown in Fig. 3, each trained to predict a distinct postoperative healing or cosmetic outcome metric derived from the finite element simulations (see Fig. 1). The first surrogate predicted breast surface deformation, quantified as the percentage of the cancer-affected breast surface exhibiting deformation greater than 1 mm at 4-week post-surgery (%BSD). The second surrogate predicted breast cavity contraction over the same 4-week period.
For validation, the high-fidelity dataset was randomly split into a training set (n = 150) and a testing set (n = 13). After training the surrogate models on the training set, predictions were generated for the testing set and compared to ground-truth simulation results. This procedure was repeated 20 times with different random splits, resulting in 260 unique test simulations evaluated against the surrogate predictions. Low-fidelity simulations were added until the error associated with the surrogates fell below 0.05 root mean square error (RMSE) for cavity contraction and 5% mean absolute error (MAE) for breast surface deformation. At this point, no further simulations were added, resulting in a total of 273 low-fidelity simulations.

Qualitative Assessment of Post-Surgical Cosmetic Outcomes
In addition to quantitative metrics, post-surgical cosmetic outcomes were qualitatively evaluated to assess the presence of visually noticeable breast surface deformations (i.e., divot or indentation). Five reviewers independently assessed the 163 high-fidelity patient simulations to determine whether a visible surface defect was present four weeks after surgery. Each reviewer, in a blinded fashion, compared preoperative breast surface renderings with the simulated postoperative outcomes and recorded a binary response indicating whether visible deformation was observed.
To translate the qualitative assessments into a predictive framework, reviewer responses for each patient were converted into a numerical consensus score. These scores were then used to train a GP surrogate model to predict the likelihood of visible post-surgical deformation using the same four patient-specific input features described earlier. Validation followed the same repeated split methodology used for the multi-fidelity surrogates.

Results

Patient-Specific Model Simulations Capture Cavity Healing and Breast Deformation Outcomes Following BCS
Building upon our previous work using a generalized breast geometry to simulate post-BCS healing [18], we applied the computational model to perform 273 low-fidelity simulations (generalized geometry) and 163 high-fidelity simulations based on patient-specific breast geometries. These simulations incorporated clinically relevant variations in breast density, breast volume, cavity volume, and cavity depth. Fig. 4 presents detailed finite element results for one representative patient, illustrating the multi-scale, spatiotemporal cavity healing response simulated over 4 weeks following BCS. The simulation aligns with trends observed in both our prior computational work and experimental porcine lumpectomy studies [11, 18]. In line with model calibration to porcine lumpectomy data, fibroblast and collagen densities at the cavity center were minimal immediately following surgery, consistent with seroma or hematoma formation [11, 18]. Cytokine concentrations within the cavity peaked and subsequently declined over the 4-week period, consistent with the early postoperative inflammatory response and its gradual resolution during healing [53, 54]. Fibroblast density increased sharply after week 1 and plateaued by week 4 (Fig. 4), coinciding with progressive deposition and remodeling of collagen by fibroblasts and myofibroblasts.
Cavity contraction was most pronounced during the first postoperative week, reaching a minimum volume shortly thereafter, before slightly rebounding and stabilizing by week 4 (Fig. 4), a pattern consistent with both porcine and clinical data [11, 55]. Permanent deformation () developed across the cavity region and adjacent healthy tissue, resulting in visible breast surface deformation (Fig. 4). At the time of tumor excision (t = 0 week), there was no change in tissue volume across the breast domain ( = 1). Shortly after surgery, permanent contracture () developed at the tissue–cavity interface, causing moderate breast surface deformation by week 1 (Fig. 4). By week 4, extensive permanent contracture was evident throughout the cavity domain, while surrounding tissue experiences tensional forces () oriented perpendicular to the cavity surface (Fig. 4). This differential deformation resulted in pronounced breast surface deformation adjacent to the cavity (Fig. 4). These findings are consistent with established tissue repair and scar formation mechanisms, where fibroblasts and myofibroblasts within the cavity contract and reorient newly deposited collagen fibers while exerting tensile forces on the surrounding ECM [36, 56].
A comparison of patient-specific simulation outcomes for breast surface deformation and cavity contraction, shown in Fig. 5, highlights the impact of individual patient characteristics, including breast density (Fig. 5A), cavity volume (Fig. 5B), breast volume (Fig. 5C), and cavity depth (Fig. 5D). Patients with lower breast density exhibited more pronounced surface deformation, including a more pronounced surface indentation, compared to those with higher breast density (Fig. 5A). Lower breast density was also associated with more rapid and severe cavity contraction (Fig. 5A). Larger cavities led to significantly greater surface deformation (Fig. 5B). While larger cavities experienced greater overall contraction, smaller cavities contracted more rapidly during the early healing phase (Fig. 5B). When evaluating the effect of breast volume, patients with larger breast volumes experienced greater magnitudes of surface deformation, however, the %BSD was similar across breast sizes (Fig. 5C). Cavity contraction followed similar trends in both cases, though patients with smaller breast volumes demonstrated slightly greater contraction by week 4 (Fig. 5C). The impact of cavity depth was also evident, with more superficial cavities leading to greater surface deformation compared to deeper cavities (Fig. 5D). Although both cases exhibited similar cavity contraction trends, patients with more superficial cavities showed marginally greater contraction by week 4 (Fig. 5D).
Breast surface deformation and cavity contraction data were extracted from each patient-specific simulation and used to train the GP surrogate models. A summary of the 163 patient-specific simulation outcomes is shown in Fig. 6, highlighting key metrics related to postoperative healing and cosmetic outcomes following BCS. On average, breast cavity volume decreased rapidly during the first 9-day post-surgery, contracting to approximately 25.70% of its original volume (Fig. 6A). The shaded 95% credible interval during this period ranged from 18.61% to 34.83% (Fig. 6A), indicating moderate variability across patients. Following this initial contraction phase, cavity volume steadily increased and stabilized, with a broader credible interval reflecting greater inter-patient variability (Fig. 6A). By week 4, the average cavity volume was approximately 51.64% of the originally excised volume, with the distribution of cavity contraction values at week 4 shown in Fig. 6B. For breast surface deformation, a significant portion of patients experienced minimal to no deformities, with 55 patients predicted to have less than 1% of the breast surface deform by more than 1 mm (Fig. 6C). These predictions were consistent with the qualitative assessment, where 88 patients were unanimously judged to have no visible post-surgical deformation (Fig. 6D). In contrast, 41 patient simulations showed more than 10% of the breast surface deforming over 1 mm, indicating more pronounced contour changes and poorer cosmetic outcomes (Fig. 6C). These findings aligned with the qualitative assessment, where 38 of these patients were unanimously judged to exhibit visible surface deformation (Fig. 6D).
Using the simulation data, we calibrated and validated multi-fidelity GP surrogate models. For each of 20 iterations, the high-fidelity dataset was randomly divided into 150 training and 13 testing samples. Surrogate predictions for the 260 testing samples were compared to computational model outputs. Validation showed that the %BSD surrogate achieved an MAE of approximately 3.48%, while the cavity contraction surrogate had an RMSE of 0.046. Additionally, when quantitative surrogate predictions of breast surface deformation were compared to the reviewer assessment consensus, the model showed agreement in approximately 81.23% of cases.

Patient Characteristics Influence Post-Surgical Breast Surface Deformation
After calibrating the GP surrogate models, they were used to predict how variations in patient-specific breast and cavity characteristics affect postoperative healing and cosmetic outcomes. Specifically, the four patient characteristics, breast volume, breast density, cavity volume, and cavity depth, were varied to determine their impact on breast surface deformation. By sampling patient characteristic values within the feasibility ranges reported in Table 1, 2,500 predictive simulations were performed. These were used to evaluate two outcomes: %BSD (Fig. 7(i)) and the likelihood of visible post-surgical deformation (Fig. 7(ii)), both as functions of the four patient characteristics. The corresponding prediction standard deviations, which quantify the uncertainty in breast surface deformation surrogate predictions across the four patient characteristics, are provided in Supplementary Fig. S1.
Consistent with clinical reports [23–26, 57–59], both breast density and cavity volume significantly impacted breast surface deformation. Lower breast densities (i.e., breasts consisting of primarily fatty tissue or scattered regions of fibroglandular tissue) and larger cavity volumes were associated with greater deformation risk (Fig. 7A(i), B(i)). The contour plot in Fig. 7A(i) indicates that for patients with low breast density, cavity volume and breast density have comparable influence on surface deformation. In contrast, for patients with higher fibroglandular tissue (> 50%), deformation is more strongly influenced by cavity volume. As shown in (Fig. 7B(i)), breast volume had a lesser effect, likely due to its correlation with breast density and cavity volume [26, 58–60]. Surface deformation was more pronounced for superficial cavities (i.e., those with smaller depth values), as shown in (Fig. 7C(i), D(i)), consistent with clinical observations [61]. However, overall cavity depth had a lesser influence on surface deformation compared to breast density and cavity volume.
Trends observed in the qualitative predictions of visible surface deformations (Fig. 7(ii)) largely mirrored those from the quantitative analyses, with some notable variations in how patient characteristics influenced deformation risk. The contour plots of cavity volume versus breast density (Fig. 7A(ii)) and breast volume versus cavity volume (Fig. 7B(ii)) showed similar patterns to their quantitative counterparts (Fig. 7A(i), B(i)), indicating that larger cavity volumes and lower breast densities were associated with a higher likelihood of visibly identifiable surface deformations. While the quantitative analysis suggested that breast density had a greater impact than cavity depth on surface deformation (Fig. 7C(i)), the qualitative predictions indicated a more balanced influence between the two (Fig. 7C(ii)). Specifically, lower values of both breast density and cavity depth increased the likelihood of an observable deformation. Similarly, although cavity volume was previously identified as the dominant factor in influencing breast surface deformation (Fig. 7D(i)), the qualitative analysis revealed that for smaller cavity volumes, cavity depth had a greater influence, while cavity volume became more impactful as its size increased (Fig. 7D(ii)).

Patient Characteristics Influence the Time-Dependent Contraction of the Breast Cavity
Using the same approach applied to breast surface deformation, the cavity contraction multi-fidelity surrogate model was utilized to predict how variations in patient-specific characteristics influenced cavity contraction over time. For each characteristic, 2,500 predictive simulations were generated by sampling across the previously established parameter ranges. These results are summarized in the contour plots in Fig. 8, illustrating the time-dependent effects of patient factors on cavity contraction. Supplementary Fig. S2 shows the corresponding prediction standard deviations, highlighting the uncertainty in the time-dependent cavity contraction surrogate predictions.
Breast cavities in low breast density regions contracted more rapidly during the first postoperative week and exhibited greater overall contraction severity over time compared to cavities in higher density (fibroglandular-rich) breasts (Fig. 8A). Similarly, larger cavity volumes were associated with greater cavity volume reduction, although smaller cavities contracted at faster rates (Fig. 8B), consistent with reported trends in human wound contraction [62, 63]. As with surface deformation, the influence of breast volume on cavity contraction was variable (Fig. 8C), likely due to interactions with breast density and excised tissue percentage [26, 58–60]. In contrast, cavity depth had little effect on either the rate or extent of cavity contraction over time (Fig. 8D).

Discussion
BCS provides the advantage of preserving breast tissue and esthetics, but its outcomes can vary greatly due to the complexities of wound healing and individual patient factors. Notably, BCS is associated with a significant incidence of post-surgical breast deformities, with approximately one-third of women developing indentations or breast asymmetry [20, 23, 26, 64]. Poor cosmetic outcomes can have a lasting impact on a patient’s self-esteem and quality of life, potentially affecting emotional well-being and overall recovery [4]. Variability in surgical outcomes and their impact on patient well-being have driven interest in personalized medicine, which guides decision-making for breast surgeons by tailoring treatment options to individual patient characteristics. For instance, risk calculators and decision models applied in clinical practice can provide valuable insights to guide treatment choices based on prior patient data. One example is a risk calculator developed by Jonczyk et al. (2021), which estimates a patient’s unique risk of acute complications following BCS, oncoplastic surgery, and various mastectomy procedures, taking into account patient demographics and other risk factors [65]. Similarly, Vos et al. developed an informed decision model that estimates the likelihood of achieving a good cosmetic outcome after BCS compared to mastectomy based on the patient’s tumor-to-breast volume ratio and tumor location [19, 20]. Additionally, Pukancsik et al. (2017) created a surgical decision model for BCS, oncoplastic surgery, and mastectomy, focusing on the location and extent of tissue excised [21]. Although these decision models are based on statistical correlations, they still do not fully capture the wound healing process or how individual patient characteristics specifically influence healing outcomes. To address this gap, we incorporated individual patient MRIs into a predictive breast healing mechanobiological model that was previously calibrated against experimental and clinical data [18]. By running simulations with the calibrated model over a large patient database from the Cancer Imaging Archive [28], we sought to advance the mechanistic understanding of how patient-specific factors impact outcomes following BCS.
The computational model used captures multi-scale cell and tissue mechanics, including complex couplings between mechanobiological activity, collagen deposition and remodeling, and deformation and biomechanics of the whole breast [18, 34–36]. To date, there are a limited number of models that depict the healing of deep soft tissue wounds. For example, Garbey and collaborators developed a 2D model based on a single patient geometry using cytokine signaling concepts from Javierre et al. (2009), adapted to predict time-dependent breast shape changes after BCS but without accounting for permanent tissue remodeling [15, 16, 66]. Vavourakis et al. (2016) developed a 3D finite element model, integrating soft tissue biomechanics with mathematical descriptions of inflammation and angiogenesis introduced by Sherratt and Murray (1991) [17, 67]. Their model was calibrated with MRI and optical surface scans from four patients taken before and 6 to 12 months after surgery [17, 67]. In contrast, our framework was calibrated against porcine and clinical data with an emphasis on (myo)fibroblast activity and collagen remodeling. These couplings enable the explanation of macroscale changes in tissue mechanics and elastoplastic deformation over long periods of time. Even though we do not model a detailed inflammatory cascade or angiogenesis, we do incorporate basic inflammation dynamics through our cytokine concentration field, similar to prior cutaneous wound healing models [14, 66, 68–70].
A key contribution of this work is the application of the calibrated 3D model to over 100 patient-specific geometries derived from preoperative MRI scans available in the Cancer Imaging Archive [28], allowing for the assessment of how patient-specific characteristics influence healing and cosmetic outcomes using a multi-fidelity GP surrogate modeling approach. These MRI scans were used to simulate individualized wound healing outcomes in 163 3D finite element simulations. The patient-specific simulations were computationally expensive, often requiring several days to complete for each individual. Thus, it was difficult to generalize from patient-specific geometries to general population insights. To enhance the dataset, low-fidelity simulations were conducted using a generalized breast geometry and a coarser mesh. This ultimately allowed the use of multi-fidelity GP regression to make efficient predictions. This approach capitalizes on the correlations between low- and high-fidelity data, outperforming models trained exclusively on high-fidelity data [46]. This methodology has previously been successfully employed by our group for the modeling and prediction of tissue expansion outcomes, demonstrating its efficacy in capturing complex mechanobiological processes with reduced computational effort [71, 72].
A similar modeling and surrogate-based approach was used by Zolfagharnasab et al. (2018), who simulated breast deformation following BCS using the model proposed by Vavourakis et al. (2016) across a limited set of patient characteristics [17, 73]. However, their study was constrained by a small sample size of patient MRIs, broad assumptions about patient-specific factors, a simplified wound healing model, and limited analysis regarding the influence of patient-specific characteristics [73]. In contrast, our multi-fidelity surrogate model was trained on over 400 finite element simulations that span a wide range of anatomic and surgical variations, more accurately reflecting the heterogeneity seen in the Cancer Imaging Archive MRI dataset [28]. While our earlier work relied on a generalized breast geometry derived from population averages [18], this study advances that foundation by incorporating patient-specific geometries from preoperative MRIs. This allowed for a more detailed investigation of how individual patient features affect healing trajectories and deformation outcomes. Our analysis revealed that breast density and cavity volume were more influential than breast volume in determining breast surface deformation (Fig. 7B), with cavity volume exerting a stronger effect in patients with denser breasts, and both parameters contributing equally in those with fattier breasts (Fig. 7A). These results are consistent with prior clinical studies [23–26, 57–61], while also expanding on existing knowledge by revealing how multiple patient-specific features interact to influence post-surgical outcomes.
For our analysis, breast surface deformation was evaluated using two complementary metrics: a quantitative measure, represented through %BSD, and a qualitative assessment, evaluating the presence of visible post-surgical surface deformation. As the two metrics quantify breast surface deformation through distinct methodologies, exact numerical agreement was not expected; rather, their comparison yields insights into the same underlying outcome. The qualitative assessment offered a more intuitive and clinically relevant interpretation, with the presence of a visible divot or indentation serving as a direct and tangible indicator of surface deformation. The quantitative data (percentage of the cancer-affected breast surface exhibiting deformation greater than 1 mm at 4-week post-surgery), is an objective measure of deformation. Despite the difference in reported values for the quantitative and qualitative methodologies, the two approaches correlate well, which indicates that they capture, to some extent, the same underlying outcome from different perspectives. Among the key findings of this comparative analysis, the quantitative measure revealed that breast density and cavity volume had a greater influence on breast surface deformation than cavity depth (Fig. 7C(i) and D(i)). In contrast, the qualitative assessment showed that cavity depth was more influential than breast density in predicting visible surface deformations (Fig. 7C(ii)). Additionally, cavity depth had a greater impact than cavity volume for smaller cavities, while their influence on surface deformation was similar for larger cavities (Fig. 7D(ii)). These findings highlight the value of integrating both quantitative and qualitative metrics to provide a more comprehensive understanding of breast surface deformation, bridging the gap between computational analysis and clinical evaluation.
This study is not without limitations. While the computational model simulations incorporated patient-specific data from preoperative MRIs to characterize breast and tumor features, longitudinal or postoperative imaging was not available to corroborate model predictive simulations against individual healing outcomes. Consequently, we are limited to compare our results against previously published clinical data focusing on the effect of patient features (breast density, tumor and breast volume, and tumor location) on healing outcomes, rather than patient-specific longitudinal data. Nevertheless, the model was previously calibrated using a combination of porcine experimental data, clinical data, and in vitro data, and model parameters were also compared against various other wound healing modeling efforts [18]. In the current implementation, many parameters remain informed by previous experimental studies or literature values rather than direct patient measurements, which limits the ability to simulate individual healing dynamics. Additionally, the biologic processes represented in the model are simplified, and future extensions could incorporate mechanisms such as angiogenesis and immune responses with additional cytokine signaling to better capture the complexity of wound healing and enhance predictive power. Several challenges also emerge when considering translation of the framework into a clinical setting. For instance, our analysis focused on breast surface deformation and cavity contraction as indicators of healing and cosmetic outcomes. Expanding this framework to include additional clinically relevant metrics, such as nipple-areolar complex asymmetry, will likely provide a more comprehensive assessment of postoperative esthetics. Important patient risk factors, such as age, diabetes, and smoking status, can also have a significant impact on healing outcomes but are not currently captured in the model. Furthermore, the model assumes homogeneous breast tissue, despite its inherently heterogeneous composition, which may limit the accuracy of predicted contracture and cosmetic outcomes. The current model also does not account for radiation therapy, which is commonly administered in the weeks following BCS [74, 75], and has been associated with fibrosis and poorer cosmetic outcomes [26, 27, 76, 77]. Future extensions of the model could incorporate additional biologic processes, including angiogenesis and immune cell–cytokine interactions, further enhancing its predictive power. Ultimately, this framework may support the development and in silico testing of therapeutic approaches such as regenerative breast tissue fillers.
In conclusion, this work presents a computational mechanobiological model of cavity healing following BCS that incorporates patient-specific data to predict individual healing outcomes and evaluate the influence of key patient characteristics on post-surgical healing and cosmetic results. By integrating preoperative MRIs with finite element simulations, we assessed how variations in breast density, cavity volume, breast volume, and cavity depth affect healing dynamics and postoperative deformation. Simulation results were used to train single- and multi-fidelity GP surrogate models, enabling efficient prediction of cavity contraction and breast surface deformation across a diverse patient population. The framework advances current approaches by not only forecasting individualized healing trajectories four weeks following surgery but also capturing tissue remodeling and deformation, offering a mechanistic understanding of post-BCS outcomes. This work also lays the foundation for a mechanobiology-driven BCS decision-making tool that is envisioned to incorporate anatomic, oncologic, surgical, and patient risk factors to simulate individualized breast deformation and healing. To enable clinical application, this framework can be implemented in a pipeline that translates preoperative imaging and patient-specific factors into predictive simulations of breast healing and cosmetic outcomes. Delivered through an interactive platform where surgeons can adjust resection parameters and visualize time-dependent cavity contraction, surface deformation, and breast asymmetry, the model would support personalized surgical planning, shared decision-making, and enhance quality of life for breast cancer survivors.

Supplementary Information
Below is the link to the electronic supplementary material.