Skip to main content

ORIGINAL RESEARCH article

Front. Physiol., 30 October 2024
Sec. Computational Physiology and Medicine
This article is part of the Research Topic Multiscale cancer modeling, in silico oncology, in silico psycho-oncology and digital (virtual) twins in the cancer domain View all 3 articles

Patient-specific prostate tumour growth simulation: a first step towards the digital twin

  • 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 (ve) and the intravascular space (vp). The exchange of substances between these compartments is modeled through the KTrans variable, which represents the extravasation rate and depends on blood flow, vascular surface area, and vascular permeability, thereby providing an overview of vascularization (Langer et al., 2010) that is closely related to oxygen levels and nutrients. Therefore, thanks to these images and different segmentation techniques, it is possible to create digital replicas of the desired organs that incorporate the vascular and cellular characteristics of each individual patient (Chase et al., 2018). This enables the potential development of patient-specific models for clinical uses, incorporating their unique parameters, which significantly enhances the comprehension of the intricate and diverse nature of these diseases (Hassanzadeh et al., 2017; Litwin and Tan, 2017; Epstein et al., 2016; Viceconti, 2015; Chase et al., 2018).

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 KTrans, as a parameter that defines vascularization characteristics. This combined and integrative approach aims to enhance the ability to predict tumour growth, understand its interactions within the prostate, and forecast PSA evolution under AS conditions. A multispecies model of partial differential reaction-advection-diffusion equations coupled with the mechanics of continuous media is here presented. Additionally, we included the dynamics of PSA with cancer development. This may help to establish a connection between patient-specific PSA levels and the aggressiveness of the cancer. In this way, the growth of the tumour can be predicted and therefore help to determine the best possible stage to apply an alternative treatment in case it is necessary. This work represents a preliminary step of a PCa model that in a future may contribute to the creation of a digital twin for PCa. The approach presented here derived from a previous tumour growth model created for neuroblastoma developed by Hervas-Raluy et al. (2024).

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 KTrans, which measures the combination of bloodflow, vessel permeability and vessel surface on each voxel, i.e., it provides useful insight into the vascularisation of the prostate (Sainz-DeMena et al., 2024).

Figure 1
www.frontiersin.org

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 KTrans maps from DCE sequences. The tumour is manually segmented by radiologists. The ADC provides an estimate of cellularity, while KTrans offers insights into the vascularization of the prostate. The Python library im2mesh (Sainz-DeMena et al., 2022) is employed to reconstruct the 3D geometry of the prostate and generate the volumetric FE mesh. Subsequently, cellularity and KTrans values are interpolated to the integration point of each element in the mesh to obtain the initial parameters of the model. Additionally, a mask is derived from nosological segmentations to identify mesh elements containing tumour cells. The Gleason Score, obtained from biopsies, provides information on the degree of cell differentiation and the overall percentage of tumour cells in each cylinder. With this data, the model is initialized for computational simulations, which are subsequently validated with clinical data.The model considerate multiple constituents (healthy, tumoural cells and stroma) and the transport of oxygen, thus the growth or shrinkage is given by proliferation or death of healthy and tumoral cells and production or remodellation of stroma. All this processes are defined by reaction-advection-diffusion equations.

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):

cellularity=ADCwADCxADCwADCmin(1)

where ADCw is the ADC value of free water molecules [ADCw=3103mm2/s (Atuegwu et al., 2013)] and ADCmin is the minimun ADC value captured. In this study, to ensure consistency of cellularity values between different patients and different follow-up MRIs, the ADCmin is standardized to 0.

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 KTrans are automatically interpolated into the integration point of each element of the FE mesh to obtain the initial parameters of the model. Through lesion segmentation, a binary mask is created, corresponding to the number of elements within the mesh. Each element is assigned a value of 1 if located within the tumour region, and 0 otherwise. This procedure is termed interpolation, as it entails the transference of data from the lesion segmentation onto the finite element mesh representing the prostate (Figure 1).

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 (pt). For stroma density, it is assumed that the fraction not occupied by cells is filled by stroma. This cellularity is physiologically unfeasible to reach a value of 1, as such a measure would imply the complete absence of stroma, which is not possible for the cells to survive. Subsequently, to calculate cell densities, it is assumed that the carrying capacity (ρci) represents the maximal cell population within a fully saturated element (Equation 2). This capacity varies between healthy and tumoural tissues, with the latter assumed to possess a greater capacity due to the cells’ propensity to deform and fill luminal spaces. The carrying capacity was determined following the methodology outlined by Atuegwu et al. (2013), assuming the number of cells that would fit in an element of 1 mm3 volume (Table 1).

ρt=cellularityptρctρh=cellularity1ptρchρs=1cellularityρcs(2)

where ρt, ρh and ρs are the population densities for tumoral cells, healthy cells and stroma.

Table 1
www.frontiersin.org

Table 1. Parameters of the model.

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 (ρi). The dynamics governing the evolution of population densities are governed by the principles of mass conservation, given by advection-reaction-diffusion (Equation 3):

ρitDρi+ρiut=ρikg11+βikρciρikρciHρoθpiρikdHθdiρo(3)

where i are the three constituents (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 kg, which is the proliferation rate, ρci, which is the tissue carrying capacity, k, which is the corrected carrying capacity, that acts as a constraint to prevent cells unlimited cell proliferation, and βi, which is the cell dependent parameter being 1 for healthy cells and 0 for tumoral cells. kg is uniformly applied across all constituents, based on the premise that βi plays a crucial role in modulating growth by diminishing the proliferation rate within healthy tissue. The kd constant defines the process of death due to hypoxia. This parameter can be defined according to the aggressiveness of the cancer (Lorenzo et al., 2016). As this work is focused on low to intermediate risk patients, the chosen parameter is selected in Table 1. Both cell proliferation and cell death processes are given by the Heaviside function H, depending on the oxygen concentration ρo, where θpi and θdi represent the oxygen concentration threshold of proliferation and death, respectively, of each constituent i. Spatial saturation is also limited by the k parameter. This parameter ensures that the concentration of the constituents never exceeds 1. It also preserves necessary space for stroma, preventing cellular overpopulation, which is vital for cell survival. The species displacement, u, results from both tumour growth and its deformation as it interacts with the surrounding tissues and organs at position x and time t. As our model focuses on predicting the tumour growth in a confined manner, the diffusive term (D), i.e., cell migration, is underestimated in these equations for simplicity, as both healthy and tumoral cell migration is assumed to be low compared to other cell processes (Pienta and Loberg, 2005).

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 KTrans, and a negative term accounting for cell consumption. Transport equation is defined as (Equation 4).

ρotDoρo=KTransρboρoih,tAoiρokoi+ρoρiρci(4)

where KTrans is the oxygen extravasation parameter, which defines the supply rate; and ρbo is the oxygen concentration in blood. The coefficients Aoi and koi (where i=h,t) represent the oxygen consumption of the different constituents of the model. More specifcally, Aoi is the maximum species consumption rate and koi is the species concentration at which the total consumption term is one-half of the total consumption term. These parameters are considered equal for both healthy and tumoral cells (Table 1), as a first simplified approach that aids in the development and analysis of the model. Although there are metabolic differences between these cell types, a single parameter is a reasonable approximation for studies primarily focused on tumour growth dynamics. Since the model’s predictions are not significantly affected by variations in this parameter, this simplification is beneficial in reducing complexity and computational cost without greatly impacting the model’s accuracy. While it is widely recognized that diffusion influences the distribution of species such as oxygen within a tumour (Greenspan, 1972; Wang et al., 2015), in well-vascularised tissue active transport through blood vessels may limit diffusion (Hervas-Raluy et al., 2024). Furthermore, it is crucial to acknowledge the considerable computational resources required to compute the diffusion process. According to our calculations, the associated costs are estimated to be around twenty times higher in cost.

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, KTrans is considered the principal physical property driving this exchange, depending on the difference in concentration between PSA in the tissue and in blood. Equation 5 and Equation 6 describe the dynamic of PSA tissue (P) and PSA serum in blood (Ps), respectively.

Pt=αhρhρch+αtρtρctKTransPPsγP(5)
dPsdt=K̃TransPPsdvγsPs(6)

where αt and αh are the production ratio of PSA by tumoral and healthy cells respectively, being αt larger than αh. γ and γs are the natural decay parameters of PSA in the tissue and in blood respectively. Serum PSA in blood is the sum of the exchange between tissue and blood in each infinitesimal part of the prostate (dv), and K̃Trans is KTrans per unit volume (K̃Trans=KTransdv).

2.3 Kinematics of growth

Let Ω be the prostate in a three-dimensional space, so that a movement φ:Ω0Ωt maps a material or reference configuration Ω0 point to a current configuration Ωt by means of Equation 7.

xX,t=φX,t(7)

The model assumes that a mechanical body consist of ith constituents that are characterised by their density and share each differential volume element. This implies that these constituents deform together, presenting the same deformation gradient tensor F. The deformation gradient tensor in the theory of continuous non-linear mechanics maps a material point X from the reference configuration (Ω0) to an spatial point x in the current configuration (Ωt) at any given time t (Equation 8).

F=xX(8)

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

F=FeiFgi(9)

where the total deformation gradient tensor F is accounted for by an elastic deformation gradient tensor (Fei), related to the stress response of the material, and an inelastic deformation gradient tensor (Fgi), connected to the volumetric growth. Therefore, to attain the current configuration Ωt from the reference configuration Ω0, an intermediate configuration related to volumetric growth is first calculated (Ωg). This intermediate configuration is not a physical state actually experienced by the tissue, it is usually an incompatible configuration (Comellas et al., 2018). The current configuration is then reached through the elastic deformation gradient component Fe, which defines the mechanical response of the tissue (Figure 2).

Figure 2
www.frontiersin.org

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 Ω0 to the current configuration Ωt involves the initial calculation of an intermediate configuration (Ωg) through Fg associated with volumetric growth. This intermediate configuration is not a physically experienced state by the tissue; rather, it represents an incompatible configuration. Subsequently, the current configuration is attained through the elastic component Fe, which characterizes the mechanical response of the tissue.

The growth component Fgi is set to be homogeneous and isotropic so that Equation 10.

Fgi=λgiI(10)

where λgi is the growth stretch ratio of every ith constituent. This factor governs inelastic deformation, which is assumed to be due to changes in mass production over the tissue corrected carrying capacity (ρci), in such a way that the growth takes place if ρgi>ρ0i and, on the contrary, resorption if ρgi<ρ0i [Equation 11 (Hervas-Raluy et al., 2024)], being ρ0i and ρgi the constituent densities in the reference and incompatible configuration, respectively.

λgi=1+kd1ρgiρiρi  if  ρgiρ0ikv1+kv21+expρgikρci  if  ρgi>ρ0i(11)

where kv1 and kv2 are constants defining the growth rate and kd1 governs the change of volume due to mass resorption.

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 Ωg to Ωt, the formulations for the densities of the grown mass (mi) concerning different configurations are expressed in Equation 12 and Equation 13.

dmi=ρgidvgi(12)
dmi=ρidv(13)

where vgi is the volume in Ωg and v the volume in Ωt. Therefore, ρgi can be rewritten as in Equation 14, where Jei is the volume ratio given by Equation 15.

ρgi=ρiJei(14)
Jei=detFei=dvdvgi(15)

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).

Ej=Eiρiρci(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 (Jgi=(1+αΔT)3, being α a coefficient expansion ratio of 1 and ΔT the increment of temperature) (Hervas-Raluy et al., 2024). By neglecting diffusion and decoupling the convective term, Equation 3 can be simplified by transforming the PDE into an ordinary differential equation (ODE), which is solved by the explicit Euler method programmed in Python. On the other hand, the Newton-Raphson method with an implicit iterative method is used to solve the growth modulus by means of Ansys®.

Figure 3
www.frontiersin.org

Figure 3. Schematic of the implementation of the model in one iteration Δt. First, the resulting oxygen concentration in each element at the beginning of the iteration is calculated with the transport equation for a Δtd (1). Since Δtd is much smaller than Δt, it is assumed that a stationary state is reached that lasts until the end of the iteration. Subsequently, the biological growth of each of the i constituents in Δt is updated taking into account the previously calculated oxygen concentration, reaching the intermediate incompatible configuration (Ωg) (2). To solve the mechanical problem (3), the software Ansys®Academic Research Mechanical, Release 19.2 is used, thus reaching the current configuration (Ωt). Finally, once the population densities in Δt are known, the PSA in tissue and blood are estimated (4). The cells in the figure represent the typical structure formed by glandular cells in the prostate, which are arranged to form lumens that secrete prostatic fluid to fulfill their glandular function. However, tumour cells proliferate, occupying and obstructing these lumens (darker cells in the figure), thereby preventing the secretion of PSA into the seminal ducts and causing a larger portion of PSA to leak into the blood. To calculate PSA dynamics, the total PSA production or decay (P*) during the time increment Δt (4.1) is first estimated. Then, the serum PSA concentration in the blood (Ps) is calculated at the stationary state (4.2). Eventually, the required tissue PSA level for equilibrium with the specified serum PSA level is recalculated (P) (4.3). When the iteration has been completed, it is checked if the final time (tfinal) of the simulation has been reached, in which case the results will be obtained. If not, a new iteration is restarted, taking as initial data the final data of the previous 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 (Kstiffness) are set in the direction normal to the prostate as contour conditions.

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 (Δt) given for each iteration of the simulation.

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 Δtd, is required for this exchange to transpire in practice. Once equilibrium is established, cellular processes are calculated and population densities are updated. The resolution of this equation employs the 8th-order explicit Runge-Kutta numerical method for ODEs. The initialisation of the oxygen concentration begins with the application of the transport equation (Equation 4), which assumes that the initial oxygen concentration in the prostate is zero. This premise facilitates the modelling of oxygen transfer from the vasculature to the tissue at the beginning of the first step.

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 Δt, resulting in a stationary state (Figure 3). To achieve this, the total amount of PSA produced or decayed (P*) in that time increment Δt is first estimated (Equation 17). Subsequently, the PSA serum concentration in the blood at the point of stationary state is calculated based on the assumed P* (Equation 18). Finally, the PSA level required in the tissue to attain the specified serum PSA level in the blood and maintain equilibrium is recalculated (Equation 19). In this way, the computation can be simplified by using the explicit Euler method for ODEs.

ΔP*Δt=αhρhρch+αtρtρctγP(17)
ΔPsΔt=KTransP*PsγsPs=0(18)
ΔPΔt=αhρhρch+αtρtρctKTransPPsγP=0(19)

Although the implementation assumes an equilibrium between blood and tissue, this does not imply that the values of P and Ps are equal at equilibrium, as there are multiple factors that influence this equilibrium. These equations are proposed on the basis of heterogeneity in KTrans, which reflects the premise that the exchange between tissue and blood is heterogeneous and vascular-dependent. The KTrans parameter summarises this idea, indicating that elevated KTrans values correspond to an enhanced exchange capacity of substances within a specific region. The initial tissue PSA is computed using Equation 19 which recalculates the tissue PSA from the serum PSA concentration. PSA serum values, derived from biochemical analyses for all patients, are adjusted via an exponential curve. The initial tissue PSA is determined by the PSA serum value on the fitted exponential curve corresponding to the date of the first MRI. This approach guarantees consistency and accurate reflection of the serum PSA measurements in the initial tissue PSA values.

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 (kg), the corrected carrying capacity (k), the incremental growth factor (kv2), the stiffness of the surrounding tissue (Kstiffness) or those related with PSA dynamics, are particularly specific to the model or display considerable variability across existing research, making difficult to assign definitive values. To address this, an optimization process is performed by means of the Python Optuna library (Akiba et al., 2019). Optuna is a widely-used hyperparameter optimization framework that leverages advanced techniques to optimize machine learning models. It primarily employs the Tree-structured Parzen Estimator (TPE) as its default sampler, utilizing a Bayesian optimization algorithm (Akiba et al., 2019). A multi-objective optimization process is chosen as it allows for adjusting the best parameters for several patients simultaneously, balancing the objective functions. Following examples from Optuna’s documentation, the optimization algorithm is developed in Python. This algorithm runs the tumour growth model a specified number of times, referring to each iteration as a trial. In each trial, the parameters to be adjusted are varied within a defined space, and the model is run for each patient with that set of parameters. Therefore, the final objective is to find the set of parameters that minimize simultaneously an objective function for each patient. Furthermore, Optuna inherently performs a sensitivity analysis throughout the optimization process, utilizing the Functional ANOVA evaluator (Akiba et al., 2019). As it adjusts parameters across trials, Optuna evaluates how variations in each parameter influence the objective function, providing insights into parameter importance.

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 (kg), the corrected carrying capacity (k), the incremental growth factor (kv2) and the stiffness of the surrounding tissue (Kstiffness). Table 2 shows the ranges within which these parameters are optimized. kg range is derived from a study in male Sprague Dawley rats where the proliferation rate of tumour cells was calculated as a function of the concentration of testosterone-activated androgen receptor (AR) (Jain et al., 2011). The k range is determined by considering the upper limit of cellularity permissible within a biological context for a given element. This ensures that the element’s cellularity remains within a physiologically realistic range. The estimation of kv2 is conducted by defining the growth threshold of an element, which is directly influenced by the density of its constituent components. This approach ensures that the element’s expansion is quantitatively aligned with its internal density parameters. Furthermore, the bladder and rectum, which are adjacent to the prostate, undergo daily fluctuations in their states. This alters their physical characteristics, posing a challenge in selecting a material stiffness that precisely mimics the surrounding anatomical conditions. The ranges given are obtained from (Dall et al., 1993; Chai et al., 2011).

Table 2
www.frontiersin.org

Table 2. Combined parameters for optimization.

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.

E=1ni=1nεoiεciεoi(20)

where i are the target values (prostate volume, tumour volume, mean prostate cellularity, and mean tumour cellularity for each follow up MRI), εoi is the objective value from clinical measurements and εci is the computational result of each trial.

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 (αt), PSA production rate of healthy cells (αh), the tissue PSA decay rate (γ) and the blood PSA decay rate (γs). Table 2 shows the ranges within which these parameters are optimized. A reduced range for PSA production rates is selected for healthy cells, predicated on the premise that in a healthy state, PSA is predominantly discharged into the prostate’s lumens and ducts, with minimal leakage into the surrounding tissue. Conversely, the proliferation of tumour cells leads to the occupation of lumen spaces, impeding the usual secretion pathways and resulting in increased PSA leakage into the tissue. Consequently, it is assumed that a small fraction of the PSA generated by healthy cells escapes into the tissue. The ranges for decay parameters γ and γs were chosen through an iterative process to optimize the model’s fit with observed PSA dynamics.

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).

MAE=1ni=1nεoiεci(21)

where i are the different time points with PSA measurements, εoi is the objective clinical value and εci is the computational result of each trial.

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
www.frontiersin.org

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
www.frontiersin.org

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 k emerges as the most influential, particularly impactful for patient B. For patient A, both k and kg exhibit significant relevance in model performance, whereas the stiffness of the environment and kv2 demonstrate comparatively lesser effects on the overall model ensemble (Supplementary Figure S1B).

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 γ emerges as the most influential, representing the decay of tissue PSA. The production of PSA of healthy cells is also significantly important (αh), followed by the production of tumoral cells (αt) and the serum PSA decay in blood (γs) (Supplementary Figure S3).

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
www.frontiersin.org

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
www.frontiersin.org

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 cm3 by the first follow up. This closely approximates the MRI observed volume of 56.36 cm3 at the corresponding time point with a relative error of 2.97% (Figure 6C). Similarly, the simulated tumour volume experiences an initial rapid increase over the first two hundred days, after which the growth rate decelerates. By the time of the first follow up, the simulated tumour volume reaches 358.36 mm3, in close agreement with the MRI observed volume of 410.15 mm3. This translates into a relative error of 15.74% (Figure 6D). Moreover, the overall cellularity in the prostate and the tumour obtained computationally is compared to that observed clinically (Figure 6E). The numerical outcome demonstrate a commendable alignment with clinical observations, with a relative error of 12.97% in the prostate and 5.66% in the tumour. The mean cellularity in the entire prostate, both clinically and computationally, hovers around 50%, while in the tumour area, it increases to approximately 60%. Even so, this correlation is rationalized by the nature of the healthy prostate, characterized as glandular tissue with numerous lumens. As tumour cells proliferate, these lumens are occupied, resulting in an elevated concentration within the tumour region. As for the PSA dynamics observed (Figure 6F), while the computational curve exhibits a slower rate of growth compared to the clinical curve, the MAE achieved is relatively low, at 0.45 ng/mL. This indicates that, despite the differences in growth rates, the computational model maintains a close approximation to the clinical data.

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
www.frontiersin.org

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
www.frontiersin.org

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 cm3. This represents a relative error of 6.15% when compared to the volume segmented from the MRI data (Figure 8C). Conversely, the tumour exhibits a notably accelerated expansion, with MRI segmentations indicating a growth rate of 80.71%. The computational model closely mirrors this progression, achieving a relative error of just 4.49% at the time of the first follow-up. Moreover, the cellularity metrics for the prostate and the tumour, as inferred from computational models, are validated against the MRI data observations (Figure 8E). The computational analysis presents a proximate reflection of the actual data, with a relative error of 6.61% for the prostate and 22.29% for the tumour. Notably, this iteration reveals a more pronounced variance between the simulated cellularity and that observed in MRI. Concerning the PSA dynamics depicted (Figure 8F), the computational model continues to project a more subdued ascent in PSA levels compared to the clinical progression, as indicated by a MAE of 1.12 ng/mL. This discrepancy underscores the model’s challenge in capturing the rapid surge in PSA levels that is evident in clinical observations.

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 KTrans, influencing the arrival of oxygen to the tissue and the exchange between PSA in the tissue and blood. Other previous prostate tumour growth models (Lorenzo et al., 2024; Jain et al., 2011) neglects this valuable patient-specific information. A limitation of the present work is that KTrans is treated as static, deviating from reality. Recognising the importance of angiogenesis in tumour growth, future efforts should explore the implementation of a dynamic KTrans for a more accurate representation. In line with earlier models (Jain et al., 2011), diffusion is neglected from PSA dynamics equations. We assume the convective transport is much more significant compared to the slower and limited impact of diffusion. The quick changes in PSA levels and the small differences in concentration further lessen the importance of diffusion (Pienta and Loberg, 2005; Tartakovsky and Dentz, 2019; Hervas-Raluy et al., 2024).

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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 [177Lu]Lu-PSMA-I&T using a population pharmacokinetic/pharmacodynamic model. EJNMMI Phys. 11, 39. doi:10.1186/s40658-024-00642-2

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

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

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu X., Reinikainen P., Kapanen M., Vierikko T., Ryymin P., Kellokumpu-Lehtinen P. L. (2018). Dynamic contrast-enhanced imaging as a prognostic tool in early diagnosis of prostate cancer: correlation with psa and clinical stage. Contrast Media Mol. Imaging 2018, 3181258. doi:10.1155/2018/3181258

PubMed Abstract | CrossRef Full Text | Google Scholar

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 Kingdom

Reviewed by:

Vasileios Vavourakis, University of Cyprus, Cyprus
Xu 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

Disclaimer: 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.