- Multiscale in Mechanical and Biological Engineering (M2BE), Aragon Institute of Engineering Research (I3A), University of Zaragoza, Zaragoza, Spain

Prostate cancer (PCa) is a major world-wide health concern. Current diagnostic methods involve Prostate-Specific Antigen (PSA) blood tests, biopsies, and Magnetic Resonance Imaging (MRI) to assess cancer aggressiveness and guide treatment decisions. MRI aligns with *in silico* medicine, as patient-specific image biomarkers can be obtained, contributing towards the development of digital twins for clinical practice. This work presents a novel framework to create a personalized PCa model by integrating clinical MRI data, such as the prostate and tumour geometry, the initial distribution of cells and the vasculature, so a full representation of the whole prostate is obtained. On top of the personalized model construction, our approach simulates and predicts temporal tumour growth in the prostate through the Finite Element Method, coupling the dynamics of tumour growth and the transport of oxygen, and incorporating cellular processes such as proliferation, differentiation, and apoptosis. In addition, our approach includes the simulation of the PSA dynamics, which allows to evaluate tumour growth through the PSA patient’s levels. To obtain the model parameters, a multi-objective optimization process is performed to adjust the best parameters for two patients simultaneously. This framework is validated by means of data from four patients with several MRI follow-ups. The diagnosis MRI allows the model creation and initialization, while subsequent MRI-based data provide additional information to validate computational predictions. The model predicts prostate and tumour volumes growth, along with serum PSA levels. This work represents a preliminary step towards the creation of digital twins for PCa patients, providing personalized insights into tumour growth.

## 1 Introduction

Prostate cancer (PCa) is a major health concern world-wide, being the most common cancer and the third leading cause of cancer-related deaths in men after lung and colorectal cancers. Projections indicate that by 2040 in Europe, the incidence rate of PCa is expected to rise by 27.6%, with mortality increasing by 53.2%, according to the World Health Organization (WHO).

The routine methods employed for the diagnosis of the PCa generally include blood tests to measure the Prostate-Specific Antigen (PSA) level, digital rectal examinations, transrectal ultrasounds, prostate biopsies, and/or imaging techniques like Magnetic Resonance Imaging (MRI) (Phan et al., 2020; Lorenzo et al., 2016; Visschere, 2018). Clinicians assign a level of aggressiveness of the cancer with these methods and establish the course of treatment to be considered (Rebello et al., 2021). Pathologists assess biopsy samples and assign a primary Gleason Score (GS), which represents the predominant histological pattern, and a secondary grade for the highest observed pattern. Both grades are assigned on a scale ranging from 1 to 5, based on microscopic architectural features and cellular characteristics. Clinicians have traditionally categorized PCa diagnoses into low, intermediate, and high-risk categories, considering a combination of Gleason patterns, PSA levels, and clinical stage (Litwin and Tan, 2017). From MRI, risk is assigned according to the PiRADs v2 protocol based on the textures of the images, the location of the tumour and its volume (Andrés et al., 2017; Weinreb et al., 2016). Specialized image acquisitions techniques can be incorporated in order to help assessing the cancer. Multiparametric MRI (mpMRI) typically includes diffusion-weighted (DWI) and dynamic contrast-enhanced imaging (DCE), in addition to T2-weighted imaging (T2w) (Litwin and Tan, 2017). Upon confirming the presence of a tumour, clinicians select the most appropriate treatment based on the patient’s risk group. Common treatments for PCa include radical prostatectomy, radiotherapy (RT), hormone therapy (HT), and active surveillance (AS). The PSA blood test plays a crucial role in various stages of PCa management, including screening, assessing future risk, detecting recurrent disease after local therapy, and managing advanced disease (Pezaro et al., 2014). AS is activated only in low-risk and some intermediate risk patients.

The integration of advanced medical imaging techniques like MRI is in line with the principles of *in silico* medicine (Rebello et al., 2021), which utilizes personalized digital models to improve disease prevention, diagnosis, prognosis, and treatment. In addition, it is possible to derive quantifiable parameters from these imaging techniques. The Apparent Diffusion Coefficient (ADC), derived from DWI, quantifies the diffusion of water molecules within tissue and it has been shown that it inversely correlates with the tissue cellularity (Atuegwu et al., 2013). Moreover, DCE-MRI sequences are employed to estimate tumour vascularization using pharmacokinetic models. The Standard Tofts Model (STM) can be applied to characterize tissue vascularization, considering two compartments: the extravascular extracellular space

Mathematical modeling in cancer research is a multifaceted field that involves critical decisions on the framework and scale of models, balancing biological accuracy with computational feasibility. Bull and Byrne (2022) identified six key mathematical hallmarks or decisions: single *versus* hybrid frameworks, homogeneity *versus* heterogeneity, spatially averaged *versus* spatially resolved, single-scale *versus* multi-scale, deterministic *versus* stochastic, and continuum *versus* discrete. Different mathematical models of cancer combine these hallmarks in different ways, leading to models that may include more biological complexity but are more challenging to analyze and/or parameterize (Greene et al. (2015); Agosti et al. (2018); Protopapa et al. (2018); Suzuki et al. (2021). For an accurate representation of tumour growth, it is essential to consider the surrounding tumour microenvironment. This includes various factors that may promote or inhibit growth. Modeling these influences often involves transport equations that account for substances like oxygen (Mpekris et al., 2015). Moreover, the model can incorporate various treatment modalities, from traditional methods like chemotherapy to cutting-edge approaches involving nanoparticles, to predict their impact on tumour progression (Dogra et al., 2019). Furthermore, the process of angiogenesis, where new blood vessels form to supply nutrients to growing tumours, is key to modeling substance delivery. Therefore, *in silico* models of tumour-induced angiogenesis are crucial for understanding and optimizing drug delivery. Moreover, angiogenesis models can be informed and tuned using routinely-acquired imaging data, e.g., angiography images, ultrasound, elastography data, or magnetic resonance images of the tumour anatomy (Hadjicharalambous et al., 2021). However, translating these models from research to clinical practice faces significant hurdles, such as parameter specification, mirroring *in-vivo* conditions, and the need for extensive validation data to ensure accuracy and reliability. Overcoming these challenges is essential for these models to enhance cancer treatment and patient outcomes (Hadjicharalambous et al., 2021).

Recent advances in PCa modeling have focused on integrating cutting-edge technologies to improve diagnosis and treatment strategies. Current models aim to describe the complex interactions between tumoral cells and their microenvironment (Phan et al., 2020). A key feature of many of these models is the inclusion of PSA as a primary indicator of PCa progression. For instance, studies by Lorenzo et al. (2016) and Mohammadi et al. (2021) used the phase-field method to account for the dynamics and co-existence of healthy and cancerous cells, as well as their interaction with nutrients, using partial differential equations (PDEs). These models are capable of predicting not only tumour growth but also PSA dynamics. Notably, Lorenzo et al. (2016) incorporated real prostate geometries extracted from CT images, enhancing the model’s personalized accuracy. Efforts have also been made to model treatment responses. Lorenzo and collaborators (2020) extended their previous work to predict tumour shrinkage under HT and Jain et al. (2011) studied the response of cancer cells under intermittent HT to predict treatment failure also taking into account PSA. In addition, PSA response to RT, factoring in tumour population and radiation-induced damage, was modeled by Lorenzo et al. (2019b).

While these works represent significant advances in PCa modeling, they still have the potential to benefit from more comprehensive clinical data obtained from routine practices, such as MRI, biopsies, and biochemical analyses. In a more recent study, Lorenzo et al. (2024) implemented a spatio-temporal mechanistic model informed by patient-specific data, using T2w and DWI MRI biomarkers to characterize tumour growth under AS. Although this is a robust and complex model, it could be further enhanced by incorporating patient-specific parameters from DCE MRI to evaluate prostate vascularization. Therefore, a more accurate representation of PCa can be obtained through mpMRI, being able not only to get an idea of the cellularity, but also of the vascularity of the prostate, a key feature in the development and treatment of this tumour.

In this work, a patient-specific model of prostate tumour growth is proposed, integrating individualised MRI-derived imaging biomarkers with a biomechanical computational model, using the ADC, as a measure of cellularity and the

## 2 Materials and methods

### 2.1 Patient-specific image biomarkers from MRI

#### 2.1.1 Patient-specific data

A retrospective study of PCa patients was conducted at the Hospital Universitario y Politécnico de La Fe de Valencia (from now on HULAFE) between 2015 and 2020, within the ProCanAid research project (PLEC 2021-007709). To be included in the study, the following criteria should be met: patients should be men over 18 years of age with a confirmed diagnosis through a positive biopsy, should have received treatment at the specified hospital, should have undergone at least two mpMRI, and should possess the necessary clinical information (GS, biopsy data, etc.). In this study, we focus on modeling tumour growth in patients under AS. A total of six patients (N = 6) meet these criteria and are included in the study. Detailed imaging acquisition information for each patient is provided in the Supplementary Material.

All patients included in this study have been diagnosed with PCa via biopsy or had a biopsy performed close to the time of their initial mpMRI scan. Each patient has a GS of (3 + 3) 6. Information on the tumour burden was obtained from the biopsy samples, providing the average percentage of tumour cells present in each biopsy cylinder for each patient. In addition, these patients have been clinically monitored with PSA measurements to track disease progression.

#### 2.1.2 mpMRI data pre-processing

Image pre-processing is performed using QP-Prostate®(Quibim S.L.) software. This software’s processing chain includes spatial smoothing, motion correction, and intra-series registration of the dynamic volumes of DWI and DCE, followed by inter-series registration of DWI and DCE to the T2W sequence. From the smoothed and motion-corrected DWI, ADC maps are calculated for each slice of the sequence. Similarly, the smoothed and motion-corrected DCE undergoes pharmacokinetic analysis (Jimenez-Pastor et al., 2023) (Figure 1). Specifically, the Standard Tofts Model (STM) is applied to the DCE and fitted using the least squares method (LSM). Among the parameters obtained, we focus on

**Figure 1**. Model initialization scheme towards the digital twin: Patient-specific data is collected from an mpMRI study, along with biopsy and biochemical reports. MRIs are analyzed using QP-Prostate®software to obtain prostate segmentation from T2w-sequences, ADC maps from DWI sequences, and

The ADC is a measure of the magnitude of diffusion (of water molecules) within tissue. It have been shown that it inversely correlates with the tissue cellularity so a quantification of the cellularity present in the tissue can be obtained by means of the Equation 1 (Atuegwu et al., 2013):

where

#### 2.1.3 Prostate and tumour segmentations

The prostate gland is automatically segmented in the T2w sequence with QP-Prostate®(Quibim S.L.) software. This software utilizes an artificial intelligence (AI) algorithm founded on Convolutional Neural Networks (CNNs).

Conversely, radiologists at HULAFE perform manual segmentation of the tumoral lesions, ascertaining the quantity of lesions discernible in the T2w sequence. In the patients selected for the study, only one lesion was found for each T2w image.

#### 2.1.4 Mesh generation and interpolation of the data

The Python library im2mesh (Sainz-DeMena et al., 2022) is used to reconstruct the 3D geometry of the different prostates and to generate the volumetric FE mesh needed for the later simulations. Im2mesh library takes prostate and tumour segmentations, makes an estimation of intermediate and missing slices, reconstruct the surface and generate a volumetric mesh. Then, the cellularity and the

#### 2.1.5 Integration of the data into the model

To integrate the image biomarkers into the model presented in the Section 2.2, additional calculations are needed. The ADC values provide a quantification of the cellularity in each element. In elements where both healthy cells and tumour cells coexist (i.e., where the tumour mask is 1), it is essential to distinguish the percentage of cellularity between these 2 cell types. This distinction is made using biopsy data on tumour load, which indicates the percentage of cellularity corresponding to tumour cells

where

### 2.2 Mathematical model

A multispecies PDE reaction-advection-diffusion model coupled with continuous media mechanics is presented for the simulation of PCa growth, encompassing both tumour growth dynamics and its impact on the entire prostate gland. This multispecies model aims to represent the phenomenological behaviour of PCa and its cellular processes, including proliferation, differentiation and apoptosis. Based on these processes, the model simulates the growth of the prostate and tumour geometry and the dynamics of PSA.

#### 2.2.1 Constituents

The prostate consists mainly of glandular epithelial cells (luminal and basal cells) and stroma, which are organised into lumens and ducts for the secretion of prostatic fluid. On the one hand, when these cells undergo mutations, they lose their healthy properties and their proliferative capacity increases, thus disrupting the formation of the glandular structure (Henry et al., 2018; Rebello et al., 2021). On the other hand, the stroma plays an essential role in the interaction between cells and the microenvironment, and its mechanical properties have been shown to play an important role in tumour growth and response (Niu and Xia, 2009). The environmental and mechanical conditions surrounding cells are essential in the control of cell populations. An example of such environmental factors would be oxygen and nutrients, normally supplied by blood vessels via diffusion (Vaupel and Kelleher, 2013; Mpekris et al., 2015; Kay, 1983). In this model, for simplicity, oxygen is considered as the only factor. This model considers mainly three constituents: healthy cells, tumoral cells and stroma; thus, each of these constituents is defined by its density

where *t* = tumoral cells; *h* = healthy cells; and *s* = stroma). The terms in the left side of the equation correspond to the temporal rate of change of the *i*-th population density, the diffusive term and the convective term and the ones in the right side correspond to the proliferation and the death term, respectively. The proliferation term depends on **x** and time t. As our model focuses on predicting the tumour growth in a confined manner, the diffusive term

#### 2.2.2 Oxygen transport

For the sake of simplicity, only oxygen and PSA (see Section 2.2.3) transport has been considered. The model can be expanded to integrate additional factors such as glucose or other nutrients. However, due to the complexity of modeling all mechanisms involved in tumour growth, oxygen has been selected as a representative factor for all these nutrients. Oxygen concentration is mainly determined by its exchange with blood, which depends specifically on prostate vascularisation represented by

where

#### 2.2.3 PSA simulation

The normal prostate architecture keeps PSA tightly confined, so that only a small proportion leaks into the circulatory system and a major part is delivered to the urethra (Sävblom et al., 2005). However, in the presence of cancer, this structure disrupts, as abnormal proliferation of tumoral cells occurs and obstructs the lumens and ducts through which PSA is secreted. Consequently, the PSA produced cannot reach the urethra, accumulating in the tissue, and leading to increased leakage into the bloodstream (Lilja et al., 2008). As a first approach, the model considers the elevation of PSA in the tissue as a consequence of the lumen disruption and leakage by tumoral cells. Therefore, the increase of PSA in the tissue is mainly proportional to the concentration of tumoral and healthy cells with different effects. The intravasation from PSA tissue to blood depends on the prostate vascularization. Thus,

where

### 2.3 Kinematics of growth

Let

The model assumes that a mechanical body consist of

Additionally, the multiplicative decomposition of the deformation tensor is used to describe the kinematics of growth (Equation 9).

where the total deformation gradient tensor

**Figure 2**. Schematic representation of the multiplicative decomposition applied to the deformation gradient **F** within a continuous framework to characterise growth: The transition from the reference configuration

The growth component

where *ith* constituent. This factor governs inelastic deformation, which is assumed to be due to changes in mass production over the tissue corrected carrying capacity

where

The consideration of mass change is essential for accurately establishing stress-strain relations to replicate tumour growth. As the mass need to be conserved from

where

All the constituents are assumed to be linear elastic and the properties of each element (j) are estimated by the rule of mixtures (Equation 16).

### 2.4 Implementation

#### 2.4.1 Numerical strategy

The proposed model is solved using the FE method. The mechanical analysis is evaluated decoupled from the biological one assuming growth incompatibility. Therefore, the growth tensor is programmed in Python 3.7.12, while the elastic contribution is computed in the commercial FE software Ansys®Academic Research Mechanical, Release 19.2, implemented by means of an ANSYS APDL (Figure 3). The thermoelastic expansion equations are used as an analogy governing volumetric changes in expansion and contraction processes to simulate tumour or prostate growth or shrinkage (Vujošević and Lubarda, 2002), so the growth or degrowth would be indicated in the ANSYS APDL file as a increment of temperature positive or negative, respectively (

**Figure 3**. Schematic of the implementation of the model in one iteration

It is assumed that each volume element contains a mixture of two structurally significant components, namely, the cells and the stroma, so that the interaction between cells and stroma is modelled. The cell populations are categorized into healthy and tumour cells, coexisting within the tumour, while the remainder of the prostate consists solely of healthy cells. Therefore, the model operates on the assumption that cellular and stromal-based matrix support the same mechanical strains, therefore as they have different mechanical properties they present different stresses. Consequently, the total stress supported by the prostate is considered the sum of both cellular and stromal contributions. To represent this, the prostate volume is discretised into two overlapping meshes with the same number of elements and sharing nodes. The mesh is composed of three-dimensional linear tetrahedral linear elements generated using im2mesh Pyhton library (Sainz-DeMena et al., 2022). The average size element is 2.5 mm for all the patients and a time step of 10 days is set.

The prostate is assumed to be located in the space bounded by the surrounding organs, mainly rectum and bladder. These organs undergo daily fluctuations due to their normal functioning, which alter their physical characteristics and make it challenging to select a material stiffness that precisely mimics the surrounding anatomical conditions. However, for this purpose, the material behaviour of these organs can be simplified to follow elastic behaviour (Dall et al., 1993; Chai et al., 2011). Therefore, to replicate the stiffness of these surrounding tissues, springs

#### 2.4.2 Multiscale temporal implementation

Cellular processes such as proliferation or death occur over a broader time span than other phenomena such as the transport of substances like oxygen and the leakage of PSA into the blood. Therefore, in order to fully integrate all these phenomena and the interaction between the cells and these substances, a multiscale temporal algorithm has been modelled. For both oxygen transport and PSA leakage, it is assumed that the equilibrium state of these substances is reached before the time increment

First, the arrival of oxygen is simulated until an equilibrium state is reached (Figure 3, step 1). By removing the diffusion term from Equation 4, an immediate equilibrium between concentrations would theoretically be attained. Nonetheless, a small time interval, denoted as

To address the multiscale temporal problem associated with PSA dynamics, it is also assumed that the equilibrium between PSA tissue concentration and PSA serum in blood is achieved well before the completion of the time increment

Although the implementation assumes an equilibrium between blood and tissue, this does not imply that the values of

### 2.5 Acquisition of model parameters

The parameters used in the model are summarised in Table 1. Most of these parameters are taken from the literature. Nevertheless, certain parameters, like the proliferation rate

Two optimization processes are performed to comprehensively refine the model. Initially, parameters governing PCa growth undergo optimization with the objective of minimizing the disparities between the MRI data and computational outcomes in prostate volume, tumour volume, mean prostate cellularity, and mean tumour cellularity. Following the acquisition of optimal values, a secondary optimization is conducted to determine parameters that most accurately align with PSA dynamics.

#### 2.5.1 PCa growth model parameters optimization

The parameters to optimise are the proliferation rate

MRI patient data serve as the basis for model calibration, encompassing measurements such as prostate volume, tumour volume, mean prostate cellularity, and mean tumour cellularity from MRI follow-ups. The simulation results are then compared with these data, being the objective minimize any disparities between them. To facilitate this, the relative error is chosen as the objective function for minimization, as it accommodates errors measured in different units, calculated according to Equation 20.

where

#### 2.5.2 PSA dynamics parameters optimization

After determining the optimal parameters for the PCa growth model, the subsequent stage involves optimizing the PSA dynamics. The parameters to optimise are the PSA production rate of tumoural cells

Clinical measurements obtained through biochemical analysis of serum PSA in blood are used for optimizing the PSA dynamics parameters. Given the variability of these measurements due to various patient and environmental factors, the mean absolute error (MAE) is selected as the objective function for minimization, owing to its robustness against outliers (Equation 21).

where

## 3 Results

### 3.1 Patient-specific data

This section introduces the patients who satisfy the study criteria detailed in Section 2.1.1.

Patient A (Figure 4) was diagnosed at the age of 60 years. He had no previous medical history of any other type of cancer and no close family members have had PCa. 199 days after biopsy diagnosis, he underwent an MRI scan which according to the PI-RADS report indicated a lesion in the medial posterior peripheral zone of the apex of the prostate with a grade 4. The biopsy analysis indicated a GS (3 + 3) 6 with an overall mean percentage of tumour volume of the extracted cylinders of 2%. He was not treated in the first instance and was included in AS with regular PSA blood tests. A follow up MRI was performed 710 days after diagnosis MRI, and finally, as clear growth was observed, he underwent radical prostatectomy surgery.

**Figure 4**. Patient A and B clinical history: Patient A, diagnosed at 60 with grade 4 PI-RADS, had a Gleason Score of 6 (3 + 3) and 2% tumour volume; underwent prostatectomy 1349 days post-diagnosis following a single follow-up MRI at day 909. Patient B, diagnosed at 68 with grade 5 PI-RADS, had a Gleason Score of 6 (3 + 3) and 3% tumour volume; received three MRIs on days 514, 703, and 907, and due to cancer’s progression, had a radical prostatectomy on day 989, followed by HT and RT. On the right, PSA measurements for both patients are displayed alongside their fitted exponential growth curves. Bellow the timeline, the FE mesh digital reconstructions of the prostate and tumours at diagnosis and follow-up MRIs is showed. The diagnosis MRI is utilized for the initialization of the model, while the subsequent follow up MRI serves for optimization.

Patient B (Figure 4) was diagnosed at the age of 68 years. His clinical record also does not indicate the existence of any other cancer or that any family members have had PCa. Diagnosis MRI, taken 121 days after biopsy diagnosis, shows a grade 5 PI-RADS lesion in the medial posterior peripheral zone of the mid and base of the prostate. The diagnostic biopsy showed a GS of (3 + 3) 6 and an overall mean tumour volume percentage of the extracted cylinders of 3%. This patient was under active follow up with biochemical analysis of PSA in the blood for approximately 2.5 years until he underwent radical prostatectomy surgery. During this time, he underwent three follow up MRI scans. This patient also exhibited extraprostatic invasion and infiltration into the seminal vesicles, indicative of a more aggressive form of cancer. Therefore, in order to eliminate any possible remaining traces of cancer cells, subsequent to the prostatectomy, he also underwent HT and RT.

Patient C (Figure 5) was diagnosed at the age of 65 years. Subsequent reviews of his medical history revealed no personal or family history of cancer. An MRI conducted 43 days prior to the biopsy pinpointed a PI-RADS grade 5 lesion located in both the mid and base sections of the peripheral zone, affecting posterior lateral and medial aspects. The biopsy returned a GS of 6 (3 + 3), with the tumorous tissue averaging 3.5% of the sampled cores’ volume. Initially opting not to undergo immediate treatment, Patient D was enrolled in AS, which included routine PSA monitoring. After a period of 666 days from the diagnosis biopsy, a subsequent follow-up MRI scan was performed showed significant tumour progression, leading to the decision for a radical prostatectomy at day 770.

**Figure 5**. Patient C and D clinical history: Patient C, diagnosed at 65 with grade 5 PI-RADS, had a Gleason Score of 6 (3 + 3) and 3.5% tumour volume; underwent prostatectomy 770 days post-diagnosis following a single follow-up MRI at day 666. Patient D, diagnosed at 56 with grade 5 PI-RADS, had a Gleason Score of 6 (3 + 3) and 3% tumour volume; he underwent 2 MRI, one prediagnosis at day −932 and another for diagnosis at day −323, and due to cancer’s progression, had a radical prostatectomy on day 501. On the right, PSA measurements for both patients are displayed alongside their fitted exponential growth curves. Bellow the timeline, the FE mesh digital reconstructions of the prostate and tumours at diagnosis and follow-up MRIs is showed. The first MRI taken for each patient are used for the initialization of the model, while the subsequent MRI serves for validation.

Patient D (Figure 5) was diagnosed at the age of 56 years. Subsequent reviews of his medical history revealed no personal or familial history of cancer. An MRI conducted 932 days prior to the biopsy pinpointed a PI-RADS grade 5 lesion located in the mid medial section of the peripheral zone. The biopsy returned a GS of 6 (3 + 3), with the tumorous tissue averaging 13.6% of the sampled cores’ volume. Initially opting not to undergo immediate treatment, Patient D was enrolled in AS, which included routine PSA monitoring. A follow up MRI was performed 710 days after diagnosis MRI, and finally, as clear growth was observed, he underwent radical prostatectomy surgery.

Patients E and F are presented in the Supplementary Material file. All six patients share the experience of undergoing AS as their principal treatment strategy. They each have GS of 6 (3 + 3), have been subject to two or more MRI scans, and have had their PSA levels measured in blood tests before and after these imaging procedures. PSA levels exhibit significant variability, yet their trajectory can be accurately characterized by an exponential curve (Karnes et al., 2018; Siebinga et al., 2024) (Figures 4, 5). In this work, two patients were chosen to fine-tune the model parameters for accurate tumour prediction (Patient A and B), while another four were utilized to validate the model’s outcomes (Patient C, D, E, and F). A mesh refinement study was performed and, to yield a good trade off between computational efficiency and accuracy, a mesh of 17833 linear tetrahedral elements and 3829 nodes was identified for the patient A, 20,572 elements and 4274 nodes for patient B, 9185 elements and 2079 nodes for patient C and 18,262 elements and 3900 nodes for patient D. Patient E was meshed with 8804 elements and 1997 nodes and patient F with 55715 elements 10848 nodes.

### 3.2 Optimization results

To perform the optimization presented in Section 2.5, patients A and B are selected.

For the optimization of the growth model, 60 trials are conducted, wherein parameters are adjusted iteratively for each patient to minimize the error function (Supplementary Figure S1A). The Optuna optimizer identifies the best trials based on those that achieve the lowest error. In this case, the optimizer did not find a single set of parameters that minimized the error for both patients simultaneously. Therefore the trial that gives an error of 0.0887 for patient A and 0.108 is selected (further information is provided in the Supplementary Material). Parameter

For the optimization of the parameters of PSA dynamics, 30 tests are conducted, wherein parameters are adjusted iteratively for each patient to minimize the error function (Supplementary Figure S2A). The minimal absolute error obtained is 0.85 ng/mL for patient A and 0.82 ng/mL for patient B. Parameter

The parameters optimised for the model are those listed in Table 1. The results of the optimization, achieved with this specific combination of parameters that yield a better fit, are illustrated in Figures 6, 7.

**Figure 6**. Patient A optimization results: **(A)** shows the simulated geometries of the prostate and the tumour compared to the MRI ones. In **(B)** the growth of the prostate and tumour volume are represented and compared to the MRI segmented volumes. In **(C)** and **(D)**, these volumes have been represented with a bar chart for follow up 1 in order to make a clearer comparison, for the prostate and tumour, respectively. **(E)** represents the overall cellularity in the prostate observed in MRI as opposed to computational outcome. Beside, the cellularity in the tumour area is shown, also comparing the MRI data and computational outcome. Finally, in **(F)**, the simulated serum PSA is compared to the clinical observations.

**Figure 7**. Patient B optimization results: **(A)** shows the simulated geometries of the prostate and the tumour compared to the MRI ones. In **(B)** the growth of the prostate and tumour volume are represented and compared to the MRI segmented volumes. In **(C)** and **(D)**, these volumes have been represented with a bar chart for each follow up in order to make a clearer comparison, for the prostate and tumour, respectively. **(E)** represents the overall cellularity in the prostate observed in MRI as opposed to computational outcome. Beside, the cellularity in the tumour area is shown, also comparing the MRI data and computational outcome. Finally, in **(F)**, the simulated serum PSA is compared to the clinical observations.

The outcomes of the computational analysis applied to patient A are described. These outcomes include the geometry of the prostate and tumour compared to one obtained from the MRI data (Figure 6A), as well as graphical representations (Figures 6B–F). Regarding volume growth the simulated prostate and tumour volumes are compared to the segmented from MRI (Figure 6B). The simulated prostate exhibits an initial rapid growth within the initial days, followed by a stabilization phase, ultimately attaining a volume of 55.14

Patient B underwent three follow up MRI scans prior to radical prostatectomy, which provided a valuable additional insight into the dynamics of cancer growth. As for patient A, the outcomes of the computational analysis applied to patient B are described. These outcomes include the geometry of the prostate and tumour compared to one obtained from the MRI data for each follow-up time point (Figure 7A), as well as graphical representations (Figures 7B–F). Regarding volume growth the simulated prostate and tumour volumes are compared to the ones segmented from MRI (Figure 7B). Similar to the scenario with patient A, the prostate volume for patient B experiences a more rapid growth approximately in the initial 300 days, eventually reaching a plateau by the first follow up. This growth pattern aligns with the observed MRI volumes, emphasizing substantial enlargement from the diagnostic image to follow up 1, obtaining relative errors of 2.76%, 4.09%, and 6.83% for follow up 1, 2, and 3, respectively. In contrast, the growth from follow up 1 to 3 exhibits a less pronounced trend. The evolution of tumour growth follows a sigmoidal pattern, characterized by gradual expansion in the initial phase, acceleration around day 350, and subsequent stabilization around day 800. This trend agrees with MRI data observations, reflecting an 70% increase in tumour size by follow-up 3 compared to its original dimensions, obtaining relative errors of 4.58%, 3.64%, and 5.52% for follow up 1, 2, and 3, respectively. Regarding cellularity (Figure 7E), MRI data observations reveal a relatively stable cellularity in the prostate across the clinical history, whereas the tumorous region experiences a gradual increase over time. Similarly, computational results depict a stabilized cellularity across the prostate, albeit with a lesser degree of variability, achieving relative errors of 25.01%, 25.25%, and 16.67% for follow up 1, 2, and 3, respectively. Notably, the cellularity in the computational tumour exhibits a more restrained growth compared to its MRI data counterpart. The results yielded relative errors of 13.52% for the follow up 1, 4.50% for the follow up 2, and 4.55% for the follow up 3. The dynamics of clinical and serum PSA are also depicted (Figure 7F). In this simulation, the PSA levels increased more rapidly than the exponential curve fitting the clinical data. Nevertheless, a MAE of 0.25 ng/mL was achieved.

### 3.3 Validation of the model

The results of the application of the model to the patients outlined in the previous section are here presented (Figures 8, 9).

**Figure 8**. Patient C simulation results: **(A)** shows the simulated geometries of the prostate and the tumour compared to the MRI ones. In **(B)** the growth of the prostate and tumour volume are represented and compared to the MRI segmented volumes. In **(C)** and **(D)**, these volumes have been represented with a bar chart for follow up 1 in order to make a clearer comparison, for the prostate and tumour, respectively. **(E)** represents the overall cellularity in the prostate observed in MRI as opposed to computational outcome. Beside, the cellularity in the tumour area is shown, also comparing the MRI data and computational outcome. Finally, in **(F)**, the simulated serum PSA is compared to the clinical observations.

**Figure 9**. Patient D simulation results: **(A)** shows the simulated geometries of the prostate and the tumour compared to the MRI ones. In **(B)** the growth of the prostate and tumour volume are represented and compared to the MRI segmented volumes. In **(C)** and **(D)**, these volumes have been represented with a bar chart for follow up 1 in order to make a clearer comparison, for the prostate and tumour, respectively. **(E)** represents the overall cellularity in the prostate observed in MRI as opposed to computational outcome. Beside, the cellularity in the tumour area is shown, also comparing the MRI data and computational outcome. Finally, in **(F)**, the simulated serum PSA is compared to the clinical observations.

The outcomes of the computational analysis applied to patient C are described following the same pattern as in previous subsection. These outcomes include the geometry of the prostate and tumour compared to one obtained from the MRI data (Figure 8A), as well as graphical representations (Figures 8B–F). Regarding volume growth the simulated prostate and tumour volumes are compared to the clinical ones (Figure 8B). The MRI prostate volume exhibits no significant increase, whereas the simulated volume demonstrates a slight initial growth before quickly stabilizing, resulting in a relative error of 3.39% during the first follow-up (Figure 8C). Conversely, the simulated tumour volume displays an initial decrease, followed by a rapid increase, eventually aligning with the MRI tumour volume with a relative error of 7.65% (Figure 8D). Additionally, the cellularity within the prostate and the tumour, as determined computationally, is corroborated by MRI data observations (Figure 8E). The numerical analysis yields a close approximation to the empirical data, with a relative error of 2.56% for the prostate and 11.40% for the tumour. A higher variability in the simulated cellularity is noted for the tumour compared to that detected in MRI. Regarding the observed PSA dynamics (Figure 8F), the computational model delineates a more gradual increase in PSA levels when contrasted with the clinical trajectory, evidenced by a MAE of 1.3 ng/mL. The model demonstrates limitations in replicating the swift escalation observed in PSA concentrations.

The outcomes of the computational analysis applied to patient D are described. These outcomes include the geometry of the prostate and tumour compared to one obtained from the MRI data (Figure 8A), as well as graphical representations (Figures 8B–F). Regarding volume growth the simulated prostate and tumour volumes are compared to the clinical ones (Figure 8B). In the computational simulations, an initial rapid increase in prostate volume is noted, which subsequently stabilizes. Upon reaching the date of the first MRI follow-up, the simulated volume attains 62.76

Results of Patient E and F are shown in the Supplementary Material. Similar trends on the predictions were obtained for these two additional patients.

## 4 Discussion

In this work, a novel approach is presented for creating a PCa model by integrating clinical MRI data, so a full representation of the whole prostate is obtained with its cellular and vascular characteristics and the location of the tumour. On top of the personalized model construction, our approach simulates and predicts the effect of tumour growth in the prostate through the FE method, coupling the dynamics of growth and the transport of oxygen to better simulate the cross-talk between cells and tumour microenvironment (Rao (2011); Ambrosi et al. (2019); Comellas et al. (2018); Hervas-Raluy et al. (2024)). In addition, our approach includes the simulation of the PSA dynamics, which allows to evaluate tumour growth through the PSA levels in the patient, that is a standard-of-care during AS.

The parameters of the biomechanical model presented here have been predominantly gathered from existing literature, yet refinement through optimization has been also performed to adjust the model. It is assumed that the mechanisms of cell proliferation and cell death are activated according to the oxygen concentration: adequate oxygen triggers proliferation, while its absence induces cell death. The model can be expanded to integrate additional factors such as glucose or other nutrients. However, due to the complexity of modeling all mechanisms involved in tumour growth, oxygen has been selected as a representative factor for all these nutrients. Future models could incorporate such additional elements by modifying the transport equations to examine their impact on tumour growth dynamics. Furthermore, the model incorporates a widely-recognized regulatory term that reduces the rate of proliferation as cell density approaches the tissue’s carrying capacity, ensuring that cell densities remain within biologically plausible limits. This term is essential for modulating the model’s growth dynamics, as demonstrated by the optimization analysis detailed in the Supplementary Material. The assumed carrying capacity of tumour cells surpasses that of healthy cells, accounting for the glandular structure in healthy tissue. Proliferation rates are further influenced by the level of cell differentiation, with more undifferentiated tumour cells exhibiting higher growth rates. The proliferation rate parameter significantly impacts growth dynamics, especially in patient A, where it ranks as the second most important factor according to the optimization analysis (see Supplementary Material). In this model, the diffusion term is underestimated from the equations that govern population density evolution. The focus is on predicting tumour growth in a localized context during the initial stages of cancer. Consequently, for simplicity, the assumption is that the cancer has not yet metastasized (Cieślikowski et al., 2020; Klotz, 2020). This simplification is a significant limitation, as it prevents the model from simulating tumour infiltration into surrounding tissues.

The distribution of oxygen within the tumour is determined by a mass transport model. To mitigate computational costs, two assumptions are incorporated. Firstly, the diffusive process of the species is considered negligible compared to the extravasation and consumption terms (Pienta and Loberg, 2005). And secondly, recognising the substantial difference in time scales between cellular processes and transport phenomena, both models are integrated into a multiscale temporal model.

In the context of PSA dynamics, a temporal model is also taken into account, postulating the attainment of equilibrium between PSA concentration in the tissue and the blood before the conclusion of each time increment. The computation of PSA in the tissue involves considering that in the tumour area, where lumens and ducts are obstructed due to abnormal cell proliferation, the accumulation of PSA in the tissue surpasses that in the healthy region (Lilja et al., 2008; Djavan et al., 1999). To account for this effect, the proliferation range of healthy cells is chosen to be much lower than that of tumour cells. Consequently, this results in a parameter for healthy cell production that is 120 times lower in our model’s calibration. It’s important to note that this parameter falls outside the typical biological range, and as such, should be interpreted with caution. With regard to natural decay rates, it is noteworthy that although the natural decay rates observed during optimization appear to deviate from the expected natural ranges, discrepancies were also observed in other models (Jain et al., 2011), so there is no clear consensus. This discrepancy highlights the inherent challenges in identifying parameters that fit the model seamlessly. Such difficulties are not uncommon and reflect the complex nature of accurately simulating real world phenomena. The production rates of PSA are critical to this dynamic, particularly those from healthy cells, constitute a greater portion of the total prostate volume, despite the fact that PSA leakage into the tissue is more intense in the tumour zone. The most significant parameter in PSA dynamics, however, is the tissue decay rate (see Supplementary Material). Moreover, vascularization is acknowledged as the primary facilitator of the exchange between PSA in the tissue and blood (Chung et al., 2014; Wu et al., 2018). The key parameter that conveys information about prostate vascularisation is

PSA is the main biomarker currently used by clinicians to assess the progression of cancer during AS, hence it seems likely that a rise in PSA is related to the progression of the disease. However, clinical PSA measurements exhibit substantial variability, as PSA levels can significantly fluctuate due to various factors, whether intrinsic to the patient -such as elevations due to benign prostatic hypertrophy, recent manipulation of the prostate due to massage or biopsy or prostatitis (Pezaro et al., 2014; Tzelepi 2023; Maeda-Minami et al., 2023)- or external. In fact, PSA variability is significant in PCa diagnostics, with studies linking fluctuations to increased cancer risk (Maeda-Minami et al., 2023). Public databases provide insights into PSA patterns in both non-PCa and PCa patients (Karunasinghe et al., 2022). PSA levels, influenced by various factors, vary with cancer stage and are not strict cutoffs. Higher PSA levels generally indicate advanced disease, but clinical stage and Gleason Score also stratify risk (Tzelepi, 2023). Despite these variations, the trajectory of the PSA data can be accurately characterized by an exponential curve (Karnes et al., 2018; Siebinga et al., 2024) (Figures 4, 5).

In this paper, two patients are used for parameter optimization of the model (patient A and B) and four others for demonstrating the preliminary feasibility of the model (patient C, D, E and F). The four six patients initially pursued AS as the primary treatment approach, involving multiple PSA measurements and MRI scans. Despite the presence of extraprostatic extension in patient B, the inclusion of patient B in the study was deemed necessary due to the challenges associated with obtaining new patient data. Moreover, given that the model is designed to simulate the tumour within a defined boundary, it remains applicable for its intended purpose. By integrating patient-specific data obtained from MRIs and biopsies, the model accounts for inter-patient variability. This approach avoids the need for individual patient data calibration, favoring a unified parameter set applicable across the patient spectrum. Consequently, it facilitates personalized predictions that are both efficient and tailored to the input data, thereby optimizing the predictive process. Finding suitable parameters for biological models is a complex endeavor. The intricate nature of biological systems means that small changes in parameters can lead to significant differences in outcomes, making it challenging to align models with real-world data. This difficulty is a notable barrier in the development of accurate predictive models. Despite the limited number of validated cases, the computationally results obtained closely align with the clinically observed outcomes. To gain a deeper understanding of the various mechanisms involved in tumour development, it would be essential to access additional data, which are currently unavailable. It is noteworthy to emphasize the difficulty in obtaining comprehensive datasets for analyses of this nature. These are retrospective studies of patients treated at the hospital, where prioritizing minimal and non-invasive tests is imperative, ensuring the patient’s health and comfort take precedence. In regards to cellularity, the overall prostate cellularity is stable over time. In contrast, tumour cellularity increases. This variance is rationalised by the stable cell density and structural maintenance in the healthy tissue, which contrasts with the abnormal cell proliferation in the tumour tissue. In the computational results we obtain a lower dispersion for both cases, but the mean aligns within the clinical observations. The reason for this is that the model inherently promotes uniformity in prostate cellularity, representing a discernible constraint. Advancements in the model’s design and further refinement will be necessary to better mimic cellularity outcome. It would indeed be advantageous to have the capability to replicate cellular distribution by modifying or introducing new hypotheses. This would enhance the model’s adaptability and accuracy in reflecting the complex biological processes within the prostate. Moreover, the complexities involved in accurately quantifying discrepancies between the cellularity maps from the model and actual observations are acknowledged. These challenges arise from the absence of direct correspondence between the mesh from the simulated prostate and the mesh generated from MRI in the follow ups, which remain unregistered and vary in their nodes and elements. Consequently, the focus has been maintained on the primary objective of the work: to simulate tumour growth, prostate, and PSA dynamics.

It is crucial to emphasize that the data used may harbor intrinsic errors originating from image acquisition procedures and assumptions made during data preprocessing. The volumes of both the tumour and prostate are susceptible to dimensional errors stemming from various factors. These include inaccuracies in segmentation, an insufficient number of MRI slices to capture the 3D characteristics, leading to potential inaccuracies in 3D reconstruction, smoothing procedures to rectify imperfections in the reconstruction, and the meshing of the geometry. Despite these challenges, our model serves as a valuable tool for simulating and understanding key aspects of PCa dynamics, providing a foundation for further refinement and improvement in the pursuit of enhanced accuracy and predictive capabilities.

This work represents a first contribution to the development of a future digital PCa twin. Future directions include the integration of different cancer treatments such as RT and HT into the model to study the effects of these treatments on tumour growth. The core vision for the digital twin of PCa is to predict the various treatment scenarios a newly diagnosed patient might face, such as AS, HT, RT, and others. This type of model holds the potential to significantly advance clinical practice. In the future, patient-specific digital twins for PCa could become a valuable tool, enabling clinicians to predict how the disease will progress in a patient-specific manner. This would aid in determining the most appropriate treatment plan, reducing the risks of both over- and under-treatment, which can cause considerable distress, improving that way the quality of life of the patient. The development of digital twins aligns with the principles of precision medicine, marking a paradigm shift in the approach to clinical decisions for PCa patients and paving the way for more effective and targeted interventions in cancer care.

## Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

## Ethics statement

The studies involving humans were approved by the Técnica del Comité de Ética de la Investigación con medicamentos del CEIM - Hospital Universitario Y Politécnico La Fe. The studies were conducted in accordance with the local legislation and institutional requirements. Written informed consent for participation was not required from the participants or the participants’ legal guardians/next of kin in accordance with the national legislation and institutional requirements.

## Author contributions

ÁP-B: Data curation, Investigation, Methodology, Software, Validation, Writing–original draft, Writing–review and editing. JG-A: Funding acquisition, Methodology, Resources, Supervision, Writing–review and editing. MG-B: Funding acquisition, Methodology, Resources, Supervision, Writing–review and editing. MAP: Conceptualization, Funding acquisition, Methodology, Project administration, Resources, Supervision, Writing–review and editing.

## Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. We would like to express our gratitude to our collaborators from HULAFE, and Quibim S.L. for their essential contribution in retrieving and processing the patient specific clinical data required for the development of this work. This research was funded by Next-Generation EU (ProCanAid Grant No. PLEC 2021-007709). The authors also acknowledge the support of the Aragon Government through Group T50 23R.

## Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

## Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

## Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2024.1421591/full#supplementary-material

## References

Agosti A., Cattaneo C., Giverso C., Ambrosi D., Ciarletta P. (2018). A computational framework for the personalized clinical treatment of glioblastoma multiforme. *J. Appl. Math. Mech.* 98, 2307–2327. doi:10.1002/zamm.201700294

Akiba T., Sano S., Yanase T., Ohta T., Koyama M. (2019). “Optuna: a next-generation hyperparameter optimization framework,” in *Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery and data mining* (Association for Computing Machinery), 2623–2631. doi:10.1145/3292500.3330701

Ambrosi D., Amar M. B., Cyron C. J., DeSimone A., Goriely A., Humphrey J. D., et al. (2019). Growth and remodelling of living tissues: perspectives, challenges and opportunities. *J. R. Soc. Interface* 16, 20190233. doi:10.1098/rsif.2019.0233

Andrés A. T., García M. O., Galiñanes M. S. (2017). Magnetic resonance imaging of the prostate: interpretation using the pi-rads v2. *Radiologia* 59, 128–138. doi:10.1016/j.rx.2016.10.005

Atuegwu N. C., Arlinghaus L. R., Li X., Chakravarthy A. B., Abramson V. G., Sanders M. E., et al. (2013). Parameterizing the logistic model of tumor growth by dw-mri and dce-mri data to predict treatment response and changes in breast cancer cellularity during neoadjuvant chemotherapy. *Transl. Oncol.* 6, 256–264. doi:10.1593/tlo.13130

Bharatha A., Hirose M., Hata N., Warfield S. K., Ferrant M., Zou K. H., et al. (2001). Evaluation of three-dimensional finite element-based deformable registration of pre- and intraoperative prostate imaging. *Med. Phys.* 28, 2551–2560. doi:10.1118/1.1414009

Boubaker M. B., Haboussi M., Ganghoffer J. F., Aletti P.(2009). Finite element simulation of interactions between pelvic organs: predictive model of the prostate motion in the context of radiotherapy. *J. Biomechanics* 42, 1862–1868. doi:10.1016/j.jbiomech.2009.05.022

Bull J. A., Byrne H. M. (2022). The hallmarks of mathematical oncology. *Proc. IEEE* 110, 523–540. doi:10.1109/JPROC.2021.3136715

Chai X., Herk M. V., Kamer J. B. V. D., Hulshof M. C., Remeijer P., Lotz H. T., et al. (2011). Finite element based bladder modeling for image-guided radiotherapy of bladder cancer. *Med. Phys.* 38, 142–150. doi:10.1118/1.3523624

Chase J. G., Preiser J. C., Dickson J. L., Pironet A., Chiew Y. S., Pretty C. G., et al. (2018). Next-generation, personalised, model-based critical care medicine: a state-of-the art review of *in silico* virtual patient models, methods, and cohorts, and how to validation them. *Biomed. Eng. Online* 17, 24. doi:10.1186/s12938-018-0455-y

Chung M. P., Margolis D., Mesko S., Wang J., Kupelian P., Kamrava M. (2014). Correlation of quantitative diffusion-weighted and dynamic contrast-enhanced mri parameters with prognostic factors in prostate cancer. *J. Med. Imaging Radiat. Oncol.* 58, 588–594. doi:10.1111/1754-9485.12230

Cieślikowski W. A., Budna-Tukan J., Świerczewska M., Ida A., Hrab M., Jankowiak A., et al. (2020). Circulating tumor cells as a marker of disseminated disease in patients with newly diagnosed high-risk prostate cancer. *Cancers* 12, 160. doi:10.3390/cancers12010160

Comellas E., Carriero A., Giorgi M., Pereira A., Shefelbine S. J. (2018). *Modeling the influence of mechanics on biological growth*. Elsevier Inc. doi:10.1016/B978-0-12-811718-7.00002-2

Dall F. H., s. Jørgensen C., Houe D., Gregersen H., Djurhuus J. C. (1993). Biomechanical wall properties of the human rectum. a study with impedance planimetry. *Gut* 34, 1581–1586. doi:10.1136/gut.34.11.1581

Djavan B., Zlotta A., Kratzik C., Remzi M., Seitz C., Schulman C. C., et al. (1999). Psa, psa density, psa density of transition zone, free/total psa ratio, and psa velocity for early detection of prostate cancer in men with serum psa 2.5 to 4.0 ng/ml. *Urology* 54, 517–522. doi:10.1016/s0090-4295(99)00153-3

Dogra P., Butner J. D., li Chuang Y., Caserta S., Goel S., Brinker C. J., et al. (2019). Mathematical modeling in cancer nanomedicine: a review. *Biomed. Microdevices* 21, 40. doi:10.1007/s10544-019-0380-2

Epstein J. I., Zelefsky M. J., Sjoberg D. D., Nelson J. B., Egevad L., Magi-Galluzzi C., et al. (2016). A contemporary prostate cancer grading system: a validated alternative to the Gleason score. *Eur. Urol.* 69, 428–435. doi:10.1016/j.eururo.2015.06.046

Faria E. C., Ma N., Gazi E., Gardner P., Brown M., Clarke N. W., et al. (2008). Measurement of elastic properties of prostate cancer cells using afm. *Analyst* 133, 1498–1500. doi:10.1039/b803355b

Greene J. M., Levy D., Fung K. L., Souza P. S., Gottesman M. M., Lavi O. (2015). Modeling intrinsic heterogeneity and growth of cancer cells. *J. Theor. Biol.* 367, 262–277. doi:10.1016/J.JTBI.2014.11.017

Greenspan H. P. (1972). Models for the growth of a solid tumor by diffusion. *Stud. Appl. Math.* LI. doi:10.1002/sapm1972514317

Hadjicharalambous M., Wijeratne P. A., Vavourakis V. (2021). From tumour perfusion to drug delivery and clinical translation of *in silico* cancer models. *Methods* 185, 82–93. doi:10.1016/j.ymeth.2020.02.010

Hassanzadeh E., Glazer D. I., Dunne R. M., Fennessy F. M., Harisinghani M. G., Tempany C. M. (2017). Prostate imaging reporting and data system version 2 (pi-rads v2): a pictorial review. *Abdom. Radiol.* 42, 278–289. doi:10.1007/s00261-016-0871-z

Henry G. H., Malewska A., Joseph D. B., Malladi V. S., Lee J., Torrealba J., et al. (2018). A cellular anatomy of the normal adult human prostate and prostatic urethra. *Cell. Rep.* 25, 3530–3542. doi:10.1016/J.CELREP.2018.11.086

Hervas-Raluy S., Sainz-DeMena D., Gomez-Benito M. J., García-Aznar J. M. (2024). Image-based biomarkers for engineering neuroblastoma patient-specific computational models. *Eng. Comput.* 40, 3215–3231. doi:10.1007/s00366-024-01964-6

Jain H. V., Clinton S. K., Bhinder A., Friedman A. (2011). Mathematical modeling of prostate cancer progression in response to androgen ablation therapy. *Proc. Natl. Acad. Sci. U. S. A.* 108, 19701–19706. doi:10.1073/pnas.1115750108

Jimenez-Pastor A., Lopez-Gonzalez R., Fos-Guarinos B., Garcia-Castro F., Wittenberg M., Torregrosa-Andrés A., et al. (2023). Automated prostate multi-regional segmentation in magnetic resonance using fully convolutional neural networks. *Eur. Radiol.* 33, 5087–5096. doi:10.1007/s00330-023-09410-9

Karnes R. J., MacKintosh F. R., Morrell C. H., Rawson L., Sprenkle P. C., Kattan M. W., et al. (2018). Prostate-specific antigen trends predict the probability of prostate cancer in a very large u.s. veterans affairs cohort. *Front. Oncol.* 8, 296. doi:10.3389/fonc.2018.00296

Karunasinghe N., Minas T. Z., Bao B. Y., Lee A., Wang A., Zhu S., et al. (2022). Assessment of factors associated with psa level in prostate cancer cases and controls from three geographical regions. *Sci. Rep.* 12, 55. doi:10.1038/s41598-021-04116-8

Kay G. F. (1983). The nature of conditioning nutrients for human malignant melanoma cultures. *J. Cell. S. D.* 62, 249–266. doi:10.1242/jcs.62.1.249

Klotz L. (2020). Active surveillance in intermediate-risk prostate cancer. *BJU Int.* 125, 346–354. doi:10.1111/bju.14935

Langer D. L., Kwast T. H. V. D., Evans A. J., Plotkin A., Trachtenberg J., Wilson B. C., et al. (2010). Prostate tissue composition and mr measurements: investigating the relationships between adc, t2, ktrans, ve, and corresponding histologic features. *Radiology* 255, 485–494. doi:10.1148/radiol.10091343

Lekka M., Gil D., Pogoda K., Dulińska-Litewka J., Jach R., Gostek J., et al. (2012). Cancer cell detection in tissue sections using afm. *Archives Biochem. Biophysics* 518, 151–156. doi:10.1016/j.abb.2011.12.013

Lilja H., Ulmert D., Vickers A. J. (2008). Prostate-specific antigen and prostate cancer: prediction, detection and monitoring. *Nat. Rev. Cancer* 8, 268–278. doi:10.1038/nrc2351

Litwin M. S., Tan H. J. (2017). The diagnosis and treatment of prostate cancer: a review. *JAMA - J. Am. Med. Assoc.* 317, 2532–2542. doi:10.1001/jama.2017.7248

Lorenzo G., Heiselman J. S., Liss M. A., Miga M. I., Gomez H., Yankeelov T. E., et al. (2024). A pilot study on patient-specific computational forecasting of prostate cancer growth during active surveillance using an imaging-informed biomechanistic model. *Cancer Res. Commun.* 4, 617–633. doi:10.1158/2767-9764.CRC-23-0449

Lorenzo G., Hughes T. J., Dominguez-Frojan P., Reali A., Gomez H. (2019a). Computer simulations suggest that prostate enlargement due to benign prostatic hyperplasia mechanically impedes prostate cancer growth. *Proc. Natl. Acad. Sci. U. S. A.* 116, 1152–1161. doi:10.1073/pnas.1815735116

Lorenzo G., Pérez-García V. M., Mariño A., Pérez-Romasanta L. A., Reali A., Gomez H. (2019b). Mechanistic modelling of prostate-specific antigen dynamics shows potential for personalized prediction of radiation therapy outcome. *J. R. Soc. Interface* 16, 20190195. doi:10.1098/RSIF.2019.0195

Lorenzo G., Scott M. A., Tew K., Hughes T. J. R., Zhang Y. J., Liu L., et al. (2016). Tissue-scale, personalized modeling and simulation of prostate cancer growth. *Proc. Natl. Acad. Sci.* 113, E7663-E7671–E7671. doi:10.1073/pnas.1615791113

Maeda-Minami A., Nishikawa T., Ishikawa H., Mutoh M., Akimoto K., Matsuyama Y., et al. (2023). Association of psa variability with prostate cancer development using large-scale medical information data: a retrospective cohort study. *Genes. Environ.* 45, 25. doi:10.1186/s41021-023-00280-7

McKeown S. R. (2014). Defining normoxia, physoxia and hypoxia in tumours - implications for treatment response. *Br. J. Radiology* 87, 20130676. doi:10.1259/bjr.20130676

Mohammadi V., Dehghan M., Marchi S. D. (2021). Numerical simulation of a prostate tumor growth model by the rbf-fd scheme and a semi-implicit time discretization. *J. Comput. Appl. Math.* 388, 113314. doi:10.1016/j.cam.2020.113314

Mpekris F., Angeli S., Pirentis A. P., Stylianopoulos T. (2015). Stress-mediated progression of solid tumors: effect of mechanical stress on tissue oxygenation, cancer cell proliferation, and drug delivery. *Biomechanics Model. Mechanobiol.* 14, 1391–1402. doi:10.1007/s10237-015-0682-0

Niu Y. N., Xia S. J. (2009). Stroma–epithelium crosstalk in prostate cancer. *Asian J. Androl.* 11, 28–35. doi:10.1038/AJA.2008.39

Pezaro C., Woo H. H., Davis I. D. (2014). Prostate cancer: measuring psa. *Intern. Med. J.* 44, 433–440. doi:10.1111/imj.12407

Phan T., Crook S. M., Bryce A. H., Maley C. C., Kostelich E. J., Kuang Y. (2020). Review: mathematical modeling of prostate cancer and clinical application. *Appl. Sci.* 10, 2721. doi:10.3390/app10082721

Pienta K. J., Loberg R. (2005). The “emigration, migration, and immigration” of prostate cancer. *Clin. Prostate Cancer* 4, 24–30. doi:10.3816/CGC.2005.n.008

Protopapa M., Zygogianni A., Stamatakos G. S., Antypas C., Armpilia C., Uzunoglu N. K., et al. (2018). Clinical implications of *in silico* mathematical modeling for glioblastoma: a critical review. *J. Neuro-Oncology* 136, 1–11. doi:10.1007/s11060-017-2650-2

Rao I. J. (2011). Modeling of growth and remodeling in soft biological tissues with multiple constituents. *Mech. Res. Commun.* 38, 24–28. doi:10.1016/j.mechrescom.2010.11.003

Rebello R. J., Oing C., Knudsen K. E., Loeb S., Johnson D. C., Reiter R. E., et al. (2021). Prostate cancer. *Nat. Rev. Dis. Prim.* 7, 9. doi:10.1038/s41572-020-00243-0

Sainz-DeMena D., García-Aznar J. M., Pérez M. A., Borau C. (2022). Im2mesh: a python library to reconstruct 3d meshes from scattered data and 2d segmentations, application to patient-specific neuroblastoma tumour image sequences. *Appl. Sci.* 12, 11557. doi:10.3390/app122211557

Sainz-DeMena D., Pérez M., García-Aznar J. (2024). Exploring the potential of physics-informed neural networks to extract vascularization data from dce-mri in the presence of diffusion. *Med. Eng. Phys.* 123, 104092. doi:10.1016/j.medengphy.2023.104092

Sävblom C., Malm J., Giwercman A., Nilsson J.-A., Berglund G., Lilja H. (2005). Blood levels of free-psa but not complex-psa significantly correlates to prostate release of psa in semen in young men, while blood levels of complex-psa, but not free-psa increase with age. *Prostate* 65, 66–72. doi:10.1002/pros.20254

Siebinga H., de Wit-van der Veen B. J., de Vries-Huizing D. M., Vogel W. V., Hendrikx J. J., Huitema A. D. (2024). Quantification of biochemical PSA dynamics after radioligand therapy with [^{177}Lu]Lu-PSMA-I&T using a population pharmacokinetic/pharmacodynamic model. *EJNMMI Phys.* 11, 39. doi:10.1186/s40658-024-00642-2

Sosa-Marrero C., Crevoisier R. D., Hernandez A., Fontaine P., Rioux-Leclercq N., Mathieu R., et al. (2021). Towards a reduced *in silico* model predicting biochemical recurrence after radiotherapy in prostate cancer. *IEEE Trans. Biomed. Eng.* 68, 2718–2729. doi:10.1109/TBME.2021.3052345

Suzuki T., Poignard C., Chaplain M., Quaranta V. (2021) “Methods of mathematical oncology,” in *Springer proceedings in mathematics and statistics 370*.

Tartakovsky D. M., Dentz M. (2019). Diffusion in porous media: phenomena and mechanisms. *Transp. Porous Media* 130, 105–127. doi:10.1007/s11242-019-01262-6

Tzelepi V. (2023). Prostate cancer: pathophysiology, pathology and therapy. *Cancers* 15, 281. doi:10.3390/cancers15010281

Vaupel P., Kelleher D. K. (2013). Blood flow and oxygenation status of prostate cancers. *Adv. Exp. Med. Biol.* 765, 299–305. doi:10.1007/978-1-4614-4989-8_42

Viceconti M. (2015). Biomechanics-based *in silico* medicine: the manifesto of a new science. *J. Biomechanics* 48, 193–194. doi:10.1016/j.jbiomech.2014.11.022

Visschere P. D. (2018). Improving the diagnosis of clinically significant prostate cancer with magnetic resonance imaging. *J. Belg. Soc. Radiology* 102. doi:10.5334/JBSR.1438

Vujošević L., Lubarda V. (2002). Finite-strain thermoelasticity based on multiplicative decomposition of deformation gradient. *Theor. Appl. Mech.* 28, 379–399. doi:10.2298/TAM0229379V

Wang Z., Butner J. D., Kerketta R., Cristini V., Deisboeck T. S. (2015). Simulating cancer growth with multiscale agent-based modeling. *Seminars Cancer Biol.* 30, 70–78. doi:10.1016/j.semcancer.2014.04.001

Weinreb J. C., Barentsz J. O., Choyke P. L., Cornud F., Haider M. A., Macura K. J., et al. (2016). Pi-rads prostate imaging - reporting and data system: 2015, version 2. *Eur. Urol.* 69, 16–40. doi:10.1016/j.eururo.2015.08.052

Keywords: prostate cancer, *in-silico* model, patient-specific, imaging biomarkers, computational oncology, finite element method (FEM)

Citation: Pérez-Benito Á, García-Aznar JM, Gómez-Benito MJ and Pérez MÁ (2024) Patient-specific prostate tumour growth simulation: a first step towards the digital twin. *Front. Physiol.* 15:1421591. doi: 10.3389/fphys.2024.1421591

Received: 22 April 2024; Accepted: 16 October 2024;

Published: 30 October 2024.

Edited by:

Henggui Zhang, The University of Manchester, United KingdomReviewed by:

Vasileios Vavourakis, University of Cyprus, CyprusXu Huang, Nanjing University of Science and Technology, China

Copyright © 2024 Pérez-Benito, García-Aznar, Gómez-Benito and Pérez. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: María Ángeles Pérez, angeles@unizar.es