# Tissue-scale, personalized modeling and simulation of prostate cancer growth

^{a}Departamento de Métodos Matemáticos e de Representación, Universidade da Coruña, 15071 A Coruña, Spain;^{b}Department of Civil and Environmental Engineering, Brigham Young University, Provo, UT 84602;^{c}Department of Information Technology, Brigham Young University, Provo, UT 84602;^{d}Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX 78712-1229;^{e}Department of Mechanical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213

See allHide authors and affiliations

Contributed by Thomas J. R. Hughes, September 23, 2016 (sent for review August 27, 2016; reviewed by Krishna Garikipati and Ellen Kuhl)

## Significance

We perform a tissue-scale, personalized computer simulation of prostate cancer (PCa) growth in a patient, based on prostatic anatomy extracted from medical images. To do so, we propose a mathematical model for the growth of PCa. The model includes an equation for the reference biomarker of PCa: the prostate-specific antigen (PSA). Hence, we can link the results of our model to data that urologists can easily interpret. Our model reproduces features of prostatic tumor growth observed in experiments and clinical practice. It also captures a known shift in the growth pattern of PCa, from spheroidal to fingered geometry. Our results indicate that this shape instability is a tumor response to escape starvation, hypoxia, and, eventually, necrosis.

## Abstract

Recently, mathematical modeling and simulation of diseases and their treatments have enabled the prediction of clinical outcomes and the design of optimal therapies on a personalized (i.e., patient-specific) basis. This new trend in medical research has been termed “predictive medicine.” Prostate cancer (PCa) is a major health problem and an ideal candidate to explore tissue-scale, personalized modeling of cancer growth for two main reasons: First, it is a small organ, and, second, tumor growth can be estimated by measuring serum prostate-specific antigen (PSA, a PCa biomarker in blood), which may enable in vivo validation. In this paper, we present a simple continuous model that reproduces the growth patterns of PCa. We use the phase-field method to account for the transformation of healthy cells to cancer cells and use diffusion−reaction equations to compute nutrient consumption and PSA production. To accurately and efficiently compute tumor growth, our simulations leverage isogeometric analysis (IGA). Our model is shown to reproduce a known shape instability from a spheroidal pattern to fingered growth. Results of our computations indicate that such shift is a tumor response to escape starvation, hypoxia, and, eventually, necrosis. Thus, branching enables the tumor to minimize the distance from inner cells to external nutrients, contributing to cancer survival and further development. We have also used our model to perform tissue-scale, personalized simulation of a PCa patient, based on prostatic anatomy extracted from computed tomography images. This simulation shows tumor progression similar to that seen in clinical practice.

According to the World Health Organization, prostate cancer (PCa) is the second most common cancer among men worldwide (1). The data are revealing: In 2012, there was an estimate of 1,095,000 new cases and 308,000 deaths worldwide associated with this cancer, and it will be responsible for about 180,890 new cases and 26,120 deaths in the United States in 2016 (2).

In most cases, PCa is an adenocarcinoma, a form of cancer that originates and develops in epithelial tissues with glandular organization, for example, the prostatic tissue in charge of producing certain substances in semen. Adenocarcinomas arise as a result of genetic alterations in cell nuclei, such as mutations and epigenetic changes. Additionally, the influence of several environmental factors, such as diet, may promote further modifications of the genome. Tumoral cells exhibit an aberrant and competitive behavior, characterized by overproliferation and an invasive demeanor. Consequently, they may endanger the structure and the normal operation of the tissue they belong to, eventually jeopardizing the survival of the whole organism. The same tumor may present different groups of cancerous cells, with different metabolism, local environment, and aggressiveness. Malignant cells may evolve during their life cycle, and thus belong to a new group. The structure, behavior, and evolution of a tumor depend on the specific genetic alterations that may have originated it and supported its evolution (3).

The development of a prostate adenocarcinoma requires a gradual accumulation of mutations in a number of different genes, which varies from patient to patient but is usually at least seven (4). As a result of these alterations of the genome over years, an initial moderate disorder of cell behavior evolves gradually toward an advanced cancer. As the tumor develops, it becomes more malignant and cell differentiation decreases.

## Current Medical Practice

Medical practice for PCa has been developed upon the above-mentioned genetic and biological bases as well as on the accumulated experience of physicians treating this disease (4). The current medical protocols (Fig. S1*A*) include all these sources of knowledge on PCa in the form of statistics for the probability of cancer stage and treatment success. In brief, PCa is easier to cure in its early stages, before it becomes excessively aggressive and spreads out of the prostate, but, unfortunately, this disease hardly ever produces any symptom until the tumor is either very large or has invaded other tissues. Presently, the best way to combat PCa is a combination of prevention and regular screening for early detection. Regular screening for PCa usually consists of a prostate-specific antigen (PSA) test and a digital rectal examination (DRE), which are usually performed periodically in men over age 50. The PSA test is a blood test for measuring the serum level of this prostate activity biomarker, which rises during PCa. The DRE is a physical test in which the doctor palpates the rectal wall near the prostate searching for unusual firm masses, as healthy prostates are normally compliant. If either the DRE or the PSA test is positive, the patient will be asked to undergo a biopsy, an invasive procedure performed with a needle guided by transrectal ultrasound (TRUS), to obtain an average of 8 to 12 tissue cores. If cancerous cells are found in the biopsy, the structure and organization of the aberrant cells will be analyzed by a pathologist to determine the aggressiveness of the tumor, which is measured by a heuristic histopathological indicator known as the Gleason score (4). With the results of the DRE, the biopsy, and some possible medical images, such as magnetic resonances or computed tomographies (CTs), PCa guidelines (5⇓–7) are used to establish a clinical stage, that is, an estimation of how far cancer has progressed. Taking the clinical stage, PSA level, and Gleason score together, protocols are used to define the risk group for a patient and associate a series of alternative treatments. The selection of the treatment, such as surgery (radical prostatectomy), cryoablation, thermal ablation, radiation therapy, hormonal therapy, chemotherapy, or a combination of them, takes into account age, life expectancy, other aspects of the clinical history, and the personal preferences of individual patients. The chosen treatment may be of a curative or palliative nature, the former being the standard for localized PCa and the latter the common practice for advanced PCa. In some cases, if the tumor is in its very early stages, or if the risk is judged to be very low, a curative treatment may be delayed while the patient is monitored until the appropriate moment for treatment (active surveillance). Patients who are not eligible for local curative treatment, or those with a short life expectancy, may benefit from a deferred palliative treatment, aimed at tackling specific symptoms of the disease, while they continue on regular screening (watchful waiting).

## Predictive Medicine and Mathematical Oncology

Despite the number of experimental and clinical investigations, a comprehensive theoretical model into which the abundance of data can be organized and understood is lacking. The emphasis of these investigations, as well as the standard clinical practice in oncology, is based on statistical patterns, which are not sufficiently accurate for individualized diagnosis, estimation of prognosis, treatment, and follow-up.

Predictive medicine (8⇓–10) is an interdisciplinary approach whose aims are the determination of personalized (i.e., patient-specific) disease progression and optimal treatment. Methods of predictive medicine are based on mathematical modeling and computer simulation of both the disease and treatments. This new approach complements the statistical and experiential basis of diagnosis and treatment planning, which constitute current medical practice (Fig. S1). In the particular case of cancer research, mathematical oncology (11⇓–13) is a new and promising field devoted to the development of models to simulate tumor growth and cancer treatment.

Although significant progress has been made in recent decades, accurately modeling the growth of cancerous tumors remains an outstanding challenge (13⇓–15). Early tumor growth models were aimed at reproducing general features of cancer, such as, for example, abnormal proliferation or evasion of programmed cell death. With the advancement of the field, new models that encode features of specific types of cancer were proposed. Arguably, the most widely studied cancer within the context of biological physics is glioblastoma multiforme, a malignant type of brain cancer (16⇓–18).

PCa has received relatively little attention from the mathematical oncology community (19⇓–21). We argue that PCa may be an ideal candidate to explore tissue-scale, personalized modeling of cancer growth for several reasons: First, the prostate is a small organ, typically the size of a walnut; second, the growth of the tumor in a given patient may be estimated using the serum PSA concentration; and third, some PCa patients do not receive any treatment, but have their PSA monitored, which potentially may open the door to in vivo model validation. Moreover, PCa is a major health problem, and is a paradigmatic example in which a predictive model could make a difference in clinical practice because there are currently many patients being overtreated and undertreated (22). Another reason why we believe it is timely to introduce a predictive model for PCa is the emergence of medical imaging as a diagnostic tool. Medical imaging is not widely used in PCa diagnosis due to insufficient confidence in current imaging modalities. However, multiparametric MRI is emerging as a promising tool that might be used for PCa diagnosis (23, 24). In Fig. S1*B*, we illustrate how imaging information can be integrated with our predictive model.

## Personalized (i.e., Patient-Specific) Modeling and Simulation of PCa

We propose a phenomenological model for PCa that qualitatively reproduces in vitro experiments and the actual growth patterns of PCa during the early to midlate stages of the disease, when PCa is usually diagnosed and treated. The computer simulations carried out with this model suggest that prostate tumors may escape starvation and necrosis by developing a tubular shape or branching. Although the potential of PSA as a PCa biomarker remains unclear, PCa diagnosis and monitoring in clinical practice relies heavily on PSA time evolution. Thus, we also propose an equation that models the PSA dynamics within the prostatic tissue. Under certain assumptions, the integration of this equation leads to another one governing the dynamics of serum PSA, whose solution can be directly interpreted by a urologist. In addition, our PSA model may shed light on the role of PSA as a biomarker for PCa, which is currently one of the most contentious debates in the urology community. Our modeling and simulation technology also addresses new challenges posed by the tissue-scale problem. Traditional smaller-scale modeling usually ignores the geometry of the host tissue. However, experimental evidence shows that tissue geometry determines sites of branching morphogenesis and plays a key role in tissue growth (25) and, consequently, in cancer development. Accounting for tissue geometry might be especially relevant in PCa because prostate anatomy varies significantly from one person to another, and PCa almost invariably develops in the peripheral zone of the prostate, close to the prostatic capsule. Thus, we propose the use of personalized anatomy in PCa modeling. We use techniques of computational geometry to extract the precise anatomy of the prostate from medical images, as well as the initial location and shape of the tumor for a patient. Then, we perform fully 3D computer simulations that predict the growth of the tumor, including a personalized simulation of PCa evolution based on prostatic anatomy obtained from CT images.

## Mathematical Model

Most continuous mathematical models of tumor growth developed to date consider several tumor species and the host tissue as different interacting phases that compete to obtain nutrients, proliferate, and occupy more space so as to thrive in the tumoral environment (26⇓⇓⇓–30). Different malignant cell types are associated with different traits, such as whether the cells are alive or not (viable, necrotic), the oxygen concentration in their local environment (normoxic, hypoxic), and whether they exhibit a certain phenotype or acquired behavior due to genetic alterations (e.g., proliferative, invasive, resistance to a certain treatment). The dynamics of each cellular species is usually modeled with a diffusion−reaction equation. In some cases, convection is also accounted for, as well as other types of migration, for example, chemotaxis, that is, cell movement driven by gradients of dissolved substances in the cell environment, or haptotaxis, that is, cell movement driven by gradients of substances attached to the extracellular matrix. The dynamics of the nutrients, as well as any other substance with a major role in cancer progression, is included by means of additional diffusion−reaction equations.

### Derivation of the Model.

In the model presented herein, we use the phase-field method to account for the healthy-to-tumoral transition and the coupled dynamics of both the host tissue and the cancerous cells, resulting in the diffusion−reaction model in Eq. **1**. The phase-field method was initially used in materials science (31, 32) to model phase transition and equilibrium between phases, eliminating the need to track interfaces with additional equations. It has been successfully used in other fields of science and engineering (33⇓⇓–36), including biomedical engineering and mathematical oncology (37). The key idea is to define an order parameter *ϕ*, which is a measure of cellular microstructure. This order parameter varies between 0 and 1. Lower values represent healthy tissue, and higher values represent an aberrant cell organization, typical of PCa. The level set

The growth of PCa is driven by a host of hormones, growth factors, and vital nutrients, but, for the sake of simplicity, we will assume that tumor growth depends on a nutrient *σ*, mainly composed of glucose. The model could also incorporate other substances, but their regulatory effect on PCa growth can be interpreted equivalently in terms of glucose. We will assume that convection has a negligible contribution to the transport of glucose, a fact that is consistent with experimental observations (38). This assumption leads to the diffusion−reaction model in Eq. **2**. The complete model reads as**1**, **1** account for nutrient-driven growth and apoptosis (i.e., programmed cell death), respectively. We assume that the growth rate depends linearly on the local nutrient concentration (39) with the growth rate coefficient *χ*. We further assume that apoptosis follows a linear relation with the population of tumoral cells, *A* being its rate, as this is the natural response in the prostatic tissue. In Eq. **2**, *s* is the nutrient supply. On the right-hand side of the equation, the third term represents the consumption of nutrient by the tumoral cells, and the fourth term accounts for the natural decay of nutrient.

The values for the model parameters can be obtained from numerous studies (39⇓⇓⇓⇓–44) on nutrient transport within tissues and tumor growth. We will always take ^{−1}·d^{−1} and ^{−1}. However, we consider several values for the growth rate *χ* and the apoptosis rate *A*, ranging from ^{−1}·y^{−1} to ^{−1}·y^{−1} and from ^{−1} to ^{−1}, to recreate the different morphologies of PCa growth during its first stages. To make the problem computationally tractable, we take *h* is the characteristic length scale of the computational mesh and *T* is the characteristic time scale, which we have considered to be 1 y (see ref. 36 for rationale). In accordance with other models of tumor growth, we will take *λ* to *λ*, which manifests the experimental observation that tumor dynamics is slower than nutrient dynamics.

We will introduce the nutrient supply in two different fashions to simulate the experimental results regarding the tumor morphology change. Beginning with a homogeneous supply with a value ^{−1}·d^{−1}, we will perform a series of simulations that already show the shape instability that makes an initially rounded tiny mass of aberrant cells grow with a fingered pattern invading the surrounding tissue. Afterward, we will set the previous value to be the base *c*, which takes values between ^{−1}·d^{−1} and 0.2 g·L^{−1}·d^{−1}. In this case, the nutrient supply will be *c* according to the vascular network depicted in the image, so that the nutrient supply will be maximum, ^{−1}·d^{−1}, near the blood vessels and minimum, ^{−1}·d^{−1}, at a point sufficiently far from them (Fig. S2). Both ^{−1}·d^{−1}.

We introduce an additional equation that models PSA dynamics within the prostate tissue based on a measure for PSA: the tissue PSA *p*, or serum PSA concentration leaked to the bloodstream per unit volume of prostatic tissue. Hence, tissue PSA will allow us to study the spatial distribution of the sources of serum PSA concentration within the gland. Both healthy and cancerous cells secrete PSA, but the latter generally produce PSA at a much higher rate than healthy cells (4). We will denote cancerous and healthy rates of PSA production per unit volume of prostatic tissue by *p* diffuses over the prostatic tissue and decays naturally with a rate *p* represents the concentration of serum PSA that leaks to the bloodstream per unit volume, we can define serum PSA *p* over the prostatic tissue **4** over the tissue region that we are simulating and assume free-flux boundary conditions, we obtain the following ordinary differential equation:**6**, chose the values ^{−1}·cm^{−3}·y^{−1} and ^{−1}, respectively, which correlated well, as can be seen in Fig. S3. Indeed, this value of

### Numerical Analysis and Simulations.

The nonlinearity of the set of equations that comprise the model and the complex geometry of both the prostate and the initial tumor in the personalized (i.e., patient-specific) simulations demand very efficient methods of resolution in terms of memory and time of computation. To perform the numerical simulations, we introduce algorithms based on the concept of isogeometric analysis (IGA) (52, 53), an emergent and cutting-edge technology that can be seen as a generalization of finite element analysis. In IGA, the standard piecewise polynomials used in the finite element method are replaced with richer functions used in computer graphics and computational geometry, such as, for example, B-splines (52, 53), nonuniform rational B-splines (NURBS) (52⇓–54), and T-splines (55). To speed up the calculations and create a technology that facilitates obtaining solutions in a clinically relevant time, we also advocate herein dynamic mesh adaptivity techniques based on the concept of hierarchical refinement and coarsening (56, 57), which enables the adaptation of the computational mesh dynamically throughout the simulation such that the approximation functions are richer close to the tumor surface, leading to dramatic savings of memory and computation time.

To extract geometric models of the prostate and the initial tumor, several methods have been proposed to construct solid anatomic NURBS models (58⇓–60). Because the geometry of a prostate is topologically equivalent to that of torus, we use a parametric mapping method to deform a torus NURBS model to match with the personalized prostate surface obtained from CT data, as depicted in Fig. 1. Finally, we can hierarchically refine the computational mesh that represents the anatomy of the prostate to perform precise personalized simulations.

CT imaging is not sufficiently accurate to detect localized PCa, so we resorted to the prior reported TRUS data to get the location and estimated shape of the tumor. Hence, to introduce the initial tumor in the model, we approximated its geometry as an ellipsoid. In particular, we L^{2}-projected (57) a 3D hyperbolic tangent field onto the spline space that defines the anatomy of the patient’s prostate to approximate the initial tumor with the phase field, hence obtaining an ellipsoidal cancerous mass in terms of the order parameter *ϕ*.

## Results

### Growth Patterns of Localized PCa.

The morphology of PCa may vary from spheroidal tumors to tubular or fingered structures that invade the prostatic epithelium forming branches, as demonstrated by clinical practice and in vitro experimental studies involving cultures of prostate epithelial cells (61). Large-scale 2D and 3D simulations of Eqs. **1** and **2** show that the model can predict both growth types, as depicted in Fig. 2, and hence qualitatively reproduces the associated clinical and experimental results for both morphologies in refs. 61 and 62.

It has been suggested that this change in tumor morphology is promoted by a series of chemical signals, also motivated by specific mutations or epigenetic changes. However, a conclusive explanation of the mechanisms of this shift in cancer growth remains largely unknown. We have been able to activate the shape instability using two mechanisms previously suggested in the literature: first, increasing the parameters *χ* and *A* in Eq. **1**, which represents a more aggressive cancer with larger apoptosis (40, 41), and, second, reducing the value of the nutrient supply *s*, which may be understood as a situation in which the tumor is surrounded by a harsh microenvironment (13, 62).

Spheroidal or ellipsoidal growth takes place in the very first stages of PCa, when it is not excessively malignant. We conjecture that, at these early stages, the tumor will feature low values for the growth rate *χ* and the apoptosis rate *A*. In particular, taking the aforementioned default values for all of the parameters in the model, ^{−1}·y^{−1} and ^{−1}, leads to a spheroidal pattern of growth, as can be seen in Fig. 2 *B* and *C*. Conversely, the fingered morphology corresponds to a more advanced cancer in comparison with the spheroidal pattern, still within the scope of localized PCa. This tubular growth can be reproduced by our model with higher values of these two parameters, which must also be balanced; that is, they must take similar values. To perform the simulations depicted in Fig. 2 *E* and *F*, we have taken *χ* = 600 L·g^{−1}·y^{−1} and *A* = 600 y^{−1}, which appear to yield a characteristic fingered morphology.

Interestingly, with this latter selection of parameters, we initially observe spheroidal growth. Should the tumor continue to develop with this morphology, the nutrient concentration in its central region would drop to values promoting hypoxia, starvation, and necrosis (Fig. S4). However, once the tumor reaches a volume that would lead to this inner harmful environment, the shift toward the fingered morphology takes place. With this growth pattern, the nutrient concentration within the tumor is not critical any longer (Fig. S4). In particular, Fig. S5 shows this evolution in a simple experimental scenario. Whether we consider the homogeneous or the heterogeneous definition for the nutrient supply, as long as *σ*.

If we reduce the value of the nutrient supply, that is, *χ* and *A* that would produce spheroidal growth otherwise. Furthermore, in the range of values for *χ* and *A* typical of fingered growth, reducing *s* accelerates the onset of the shape instability. We also observe that a drop in the values of *s* leads to slower tumor growth, whereas increasing the nutrient supply makes the tumor grow at a faster rate. In the case of the fingered morphology, the higher the value of *s*, the thicker the branches of the tumor (Fig. S6), as the resulting higher concentrations of the nutrient within the tissue support larger tumoral structures.

### Personalized Simulation.

The Austin Radiological Association provided us with anonymized axial and coronal CT images of a small cohort of patients with PCa, with and without contrast (institutional review board approval and informed consent were not required for this study). CT is a not a suitable technique for the detection of localized PCa, due to its reduced contrast between healthy and cancerous tissue within this gland; however, it is implemented in PCa protocols primarily to scan for metastasis or as a resource to initially analyze the prostate anatomy and the pelvic region with more detail than ultrasound.

We have chosen one of these patients to predict the evolution of his prostatic tumor with our model and assess the performance of the simulation on the basis of the current knowledge of PCa evolution. This patient showed a small cancerous lesion of 0.03 cm^{3}, previously detected by ultrasound, in the posterior part of the peripheral zone at the middle region of the prostatic gland. The volume of the prostate gland is considered constant along the simulation and has a value of 71.67 cm^{3}.

Fig. 3 represents the prediction of the growth of a prostatic adenocarcinoma according to our model at tissue scale and on a personalized basis; that is, we ran our model on the actual geometry of the prostate of the patient under consideration, after extracting it and the initial tumoral volume from the CT output files. This simulation represents tumor growth during 1 y, and, at the end, the tumor volume is 2.66 cm^{3}. Additionally, we have produced Movies S1–S3 showing this predicted tumor growth from an anterior view, a posterior view, and on an axial section at the middle gland, respectively.

Initially, the tumor grows following an ellipsoidal pattern, but it soon starts to develop preferential directions of growth, defining finger-like structures. As the tumor expands within the prostate, its geometry evolves toward a structure consisting of thickened layers of cancerous tissue with fingered protrusions along the main directions of growth. The axial section at the midgland shows an evolution from the spheroidal pattern of growth to the fingered morphology similar to that observed in large-scale 2D and 3D simulations. The tumoral area in these sections qualitatively matches some delineations of PCa made on actual histopathologic specimens (63) in terms of size and shape. As the delineations of ref. 63, Fig. 3 *D*–*F* shows axial sections with disjoint cancerous regions. In the absence of the 3D representation of both the prostate and the tumor, this situation might be classified as a multifocal tumor rather than a single tumor with fingered shape.

### Tissue and Serum PSA.

In our simulations, tissue PSA evolution is obtained from Eq. **4**, and serum PSA evolution is then obtained from Eq. **6**. The field of tissue PSA typically shows higher values in the central areas of the tumor and lower values as we move away from the lesion, as depicted in Fig. 4. Additionally, Movie S4 shows how the *p* field evolves in time on an axial section of the prostate at the midgland. The difference in the values of *p* in healthy and cancerous tissue is due to the fact that the initial tumor starts producing PSA at a higher rate from the first moment, and, as the tumor grows, higher tissue PSA values are found in the areas occupied by cancerous cells and following the preferential directions of growth. In general, tissue PSA values do not differ greatly between healthy and tumoral tissue in our 2D and 3D experimental simulations. However, in the tissue-scale, personalized simulations, the values of *p* in the tumor, around 0.70 ng·mL^{−1}·cm^{−3}, were an average of 3 times higher than in the healthy regions, around 0.24 ng·mL^{−1}·cm^{−3}, as seen in Fig. 4 and Movie S4.

Taking lower values for

We obtain the serum PSA by integrating the values of *p* over the entire prostate anatomy. We lacked the biochemical data for the patients on our cohort, so we had to estimate the baseline serum PSA according to the tumor stage and the dynamics imposed by our model. Beginning with the minimum value of 4.0 ng·mL^{−1}, beyond which a physician would recommend a TRUS-guided biopsy, we observed a steep jump toward a higher value in the first time steps, and then serum PSA followed exponential dynamics, as observed in clinical practice (4, 46, 47, 64). Therefore, to produce more realistic simulations, we restarted our model once this readjustment of serum PSA was completed, correcting the initial condition for *p* so that it corresponded to the value of serum PSA reached just after the rising branch. Hence, serum PSA readjustment was eliminated, and it followed exponential dynamics from the beginning of the simulations. In particular, we have estimated an initial serum PSA level of 17.3 ng·mL^{−1}, and, at the end of the simulation, PSA reaches 18.9 ng·mL^{−1}. Fig. 4 shows the predicted serum PSA history for our patient, taking into account all these considerations together with the time history for tumor volume. Indeed, we have observed a parallelism between both dynamics, as has been previously reported in literature (4, 46, 47).

## Discussion

Predictive medicine and mathematical oncology offer a cutting-edge approach to medical practice. These new trends may be the key to organizing and broadening our understanding of some critical diseases nowadays, such as PCa, by providing robust, comprehensive theoretical models integrating all of the major known and forthcoming results from research in biological, biophysical, medical, and biomedical sciences as well as engineering and computational mechanics. In particular, the simple model presented herein is able to reproduce the typical features of prostatic adenocarcinoma growth during the first to midlate stages of this disease as seen in experimental results and clinical practice. We believe this might be the first step toward the derivation of a complete model that, in the context of Fig. S1*B*, could be implemented in clinical software for physicians to aid in the diagnosis of the disease, the prediction of its evolution, assessment of alternative treatments, and management of follow-up.

### Tissue-Scale, Personalized Simulations.

We have carried out personalized simulations (Fig. 3) over the actual anatomy of a patient, introducing the tumor morphology at the time of diagnosis by extracting both geometries from CT output data. We observed that, if we changed the position of the tumor within the prostate but conserved its original shape, the tumoral pattern of growth was different. These results suggest that cancer growth morphology is dependent on the particular geometry of the tumor as well as on the specific anatomy of the prostate. Consequently, this observation supports the major importance of considering personalized anatomies for both the prostate and the tumor in order to accurately predict tumor growth. Introducing the particular anatomy of a patient in the model is crucial for predicting the potential risks that PCa growth may entail, such as prostatic capsule penetration, invasion of nearby tissues, or tumoral expansion to the different parts of the prostate, hence triggering symptoms such as impotence or irregularities in urination and ejaculation (4). The possibility of predicting when these risks or symptoms will appear would be a major advance in watchful waiting and active surveillance. Hence, the optimization of these two clinical options to manage the disease, including the modeling and simulation of PCa, would decrease the overtreatment of patients, in many instances leaving them relieved or cured after the delivery of the corresponding treatment, with minimum side effects and a good quality of life. The inclusion of the individual’s prostate and tumor geometries is essential to adequately follow up the progression of the disease in terms of tumor volume and location and PSA production, both before and after the delivery of a certain treatment. Predicting the actual tumor shape and volume within the anatomy of the prostate of a given patient would enable optimization of surgery and radiation therapy planning and delivery, obtaining high precision and efficiency rates.

The simulation presented in Fig. 3 shows an important and typical feature of PCa, namely, why some biopsies are falsely negative or apparently yield a lower tumoral stage. As the tumor grows, it develops finger-like structures after the morphology shift that, in a 3D, personalized scenario, evolve not only as branches but also as thickened layers of cancerous cells. Likewise, there are other layers and regions of healthy tissue surrounding and in between these cancerous layers, which are required to provide an adequate distribution of the nutrient. When a biopsy is performed following the current standards, a number of tissue samples are obtained in prescribed positions trying to map the whole prostatic gland. However, the tumor might be out of the reach of the needle used in the biopsy, resulting in a false negative, or it may target one of these tumoral thickened layers orthogonally, what may produce a low percentage of tumor in the tissue sample, leading to further biopsies or/and an erroneous stage of PCa that may compromise the patient quality of life and survival. One might conclude that the pathology analysis of a prostatic tumor biopsy would always provide a lower bound to the severity of the cancer, creating a dangerous predictive scenario.

Ideally, the combination of multiparametric MRI and computational modeling and simulation of PCa may eliminate the biopsy of the prostate during PCa diagnosis and staging, as proposed in Fig. S1, hence dropping an invasive procedure with questionable outcomes and under continuous debate in the medical community. In fact, the combination of T2-weighted MRI, diffusion-weighted MRI, dynamic contrast-enhanced MRI, and, optionally, magnetic resonance spectroscopic imaging has already been able to detect and stage PCa with high sensitivity and specificity (23, 24), comparable to those of current standard biopsies. As a result, multiparametric MRI is now increasingly regarded as a future alternative to biopsy or, at least, a technique to be fused to TRUS-guided biopsy to optimize this procedure (65).

### Modeling and Simulation of PSA Kinetics.

We have also included an equation to model the behavior of PSA by introducing a measure, tissue PSA or the serum PSA concentration leaked to the bloodstream per unit volume of prostate tissue (Eq. **4**). Once integrated over the computational domain that represents the anatomy of the prostate, serum PSA dynamics may be followed as shown in Eq. **6**. The tissue PSA allows us to analyze where PSA is being leaked to the bloodstream and the amount of this PCa biomarker discharged by each type of tissue, as depicted in Fig. 4. Furthermore, the results presented in Fig. 4 also show that our model is able to compute a realistic evolution for serum PSA, similar to those seen in the literature (47). Tissue PSA links serum PSA to tumor burden in our model, hence establishing a theoretical basis for the link between both dynamics. A corroboration of this connection arose when we estimated the baseline serum PSA for our simulation. For a given selection of parameters in Eqs. **1** and **4**, and the initial tumoral volume, the initial serum PSA was intrinsically determined. Therefore, if we began the simulation with a different serum PSA value, a steep branch would appear until it was readjusted. Additionally, it can be seen in Fig. 4 that the evolution of the tumoral volume and that of serum PSA follow similar dynamics. This behavior is typically seen in low- and intermediate-risk tumors (4, 47). Nonetheless, the connection between tumor growth and serum PSA dynamics is an ongoing major debate in the urologic community. Regular screening of PCa relies on PSA tests, but, even though most prostatic adenocarcinomas grow showing a parallelism between the tumor volume and PSA levels in blood, it is known that there exists a certain degree of divergence between the dynamics of these two variables; this may even reach the limit of patients with low PSA and high tumoral masses in their prostates, or vice versa. Advanced and metastatic PCa may also exhibit almost no connection between serum PSA and tumor burden. However, PSA screening has greatly improved the early detection and treatment of PCa, leading to higher rates of incidence and an improved efficiency in its management worldwide. Therefore, our opinion is that this biomarker should be implemented in a mathematical model of PCa, as it permits direct connection with existing clinical data.

### Future Work.

There is still much room for improvement in PCa modeling and simulation. First of all, several cancerous species could be considered, for instance, normoxic, hypoxic, and necrotic, so as to better reproduce the shape instability. The model could also account for other substances that are known to have a role in PCa, such as oxygen or testosterone. Haptotaxis could play a role in more aggressive PCa, and the invasiveness of PCa could be better reproduced if the diffusive term in Eq. **1** is defined in a more sophisticated fashion. Another aspect of interest is the study of mechanical stress on tumor growth, which has been studied experimentally (66, 67) and computationally (68). This could be particularly important in PCa, because the prostate is subjected to significant mechanical confinement. PSA modeling in the medical community may be worth revisiting to find further connections between tumor shape and volume and PSA production, hence improving the definition of the rates **4**. The selection of the parameters in Eq. **4** could also be personalized by using data of a particular patient. Men over age 50 have their PSA monitored routinely. Methods of inversion, Bayesian statistics, or machine learning could be used to determine the parameters from these data. We would also want to explore the modeling and simulation of the effects of radiation therapy in PCa with our model, following previous studies carried out for low-grade gliomas (69) and glioblastoma multiforme (70). To perform full predictive tissue-scale, personalized simulations of PCa growth and assess their performance, it is mandatory to have some basic data on each patient over a large period, namely, high-quality medical images, biochemical data, and clinical reports.

## Conclusion

The work presented herein is just an initial example of the potential that the computational modeling and simulation of cancer may entail. Introducing a comprehensive theoretical model and the technology for its simulation in medical practice for PCa would provide a personalized (i.e., patient-specific) diagnosis and treatment. It might also provide a basis for noninvasive diagnosis, avoiding biopsies, and a more precise diagnosis, thus eliminating unwarranted invasive treatments. Mathematical oncology could also rationalize and organize our current and future knowledge of PCa, guiding the future interdisciplinary research efforts in the fields of biology, medical science, biophysics, and biomedical and computational engineering, so as to enrich and further validate the model for its use in actual clinical practice.

## Acknowledgments

We acknowledge Kent Beasley, MD, and the Austin Radiological Association for providing anonymized medical images and reports of PCa patients. We also thank the Fulton Supercomputing Lab at Brigham Young University for use of their facilities to run the tissue-scale, personalized simulations. H.G., G.L., and G.V. were partially supported by the European Research Council (Contract 307201) and Xunta de Galicia (Conselleria de Cultura, Educacion e Ordenacion Universitaria).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: hughes{at}ices.utexas.edu.

Author contributions: G.L., T.J.R.H., and H.G. designed research; G.L., T.J.R.H., and H.G. performed research; G.L., M.A.S., K.T., T.J.R.H., Y.J.Z., L.L., G.V., and H.G. contributed new reagents/analytic tools; G.L., T.J.R.H., and H.G. analyzed data; and G.L., T.J.R.H., and H.G. wrote the paper.

Reviewers: K.G., University of Michigan; and E.K., Stanford University.

Conflict of interest statement: T.J.R.H., H.G., and G.L. are listed as coinventors in a patent application that has been filed by The University of Texas at Austin. The invention emanates from the research in this paper.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1615791113/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵.
- International Agency for Research on Cancer, World Health Organization

*GLOBOCAN 2012: Estimated Cancer Incidence, Mortality and Prevalence Worldwide in 2012.*Version 1.2. Available at globocan.iarc.fr/Default.aspx. Accessed February 23, 2016. - ↵.
- American Cancer Society

*Cancer Facts & Figures 2016*. Available at www.cancer.org/acs/groups/content/@research/documents/document/acspc-047079.pdf. Accessed February 26, 2016. - ↵.
- Alberts B, et al.

- ↵.
- Walsh PC,
- Worthington JF

- ↵.
- National Comprehensive Cancer Network

- ↵.
- Mottet N, et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Jain HV,
- Clinton SK,
- Bhinder A,
- Friedman A

- ↵.
- Morken JD,
- Packer A,
- Everett RA,
- Nagy JD,
- Kuang Y

- ↵
- ↵
- ↵
- ↵
- ↵.
- Nelson CM,
- Vanduijn MM,
- Inman JL,
- Fletcher DA,
- Bissell MJ

- ↵.
- Bearer EL, et al.

- ↵.
- Gallaher J,
- Babu A,
- Plevritis S,
- Anderson ARA

- ↵.
- Lima EABF,
- Oden JT,
- Almeida RC

- ↵
- ↵
- ↵.
- Qin RS,
- Bhadeshia HK

- ↵.
- Qin RS,
- Bhadeshia HKDH

- ↵
- ↵.
- Gomez H,
- Cueto-Felgueroso L,
- Juanes R

- ↵.
- Bueno J,
- Bona-Casas C,
- Bazilevs Y,
- Gomez H

- ↵.
- Gomez H,
- Calo VM,
- Bazilevs Y,
- Hughes TJR

- ↵.
- Vilanova G,
- Colominas I,
- Gomez H

- ↵.
- Ciarletta P,
- Foret L,
- Ben Amar M

- ↵.
- Ellem KAO,
- Kay GF

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Gerhardt H, et al.

- ↵
- ↵.
- Vollmer RT

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Cottrell JA,
- Hughes TJR,
- Bazilevs Y

- ↵.
- Piegl L,
- Tiller W

- ↵.
- Bazilevs Y, et al.

- ↵.
- Scott MA,
- Thomas DC,
- Evans EJ

- ↵.
- Thomas DC,
- Scott MA,
- Evans JA,
- Tew K,
- Evans EJ

- ↵.
- Floater MS,
- Hormann K

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Vollmer RT,
- Humphrey PA

- ↵.
- Sandler K, et al.

- ↵
- ↵.
- Tse JM, et al.

- ↵.
- Mills KL,
- Kemkemer R,
- Rudraraju S,
- Garikipati K

- ↵.
- Pérez-García VM, et al.

- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Engineering