Computational-Model-Based Biopharmaceutics for p53 Pathway Using Modern Control Techniques for Cancer Treatment

Abstract

The p53 pathway has been the focus of many researchers in the last few decades owing to its pivotal role as a frontline cancer suppressant protein. It plays a vital role in maintaining cell cycle checkpoints and cell apoptosis in response to a broken DNA strand. This is why it is found in the mutated form in more than 50% of malignant tumors. To overcome this, various drugs have been proposed to revive the p53 pathway in cancer patients. Small-molecule-based drugs, such as Nutlin 3a, which are capable of performing this stimulation, are at the fore of advanced clinical trials. However, the calculation of their dosage is a challenge. In this work, a method to determine the dosage of Nutlin 3a is investigated. A control-systems-based model is developed to study the response of the wild-type p53 protein to this drug. The proposed strategy regulates the p53 protein along with negative and positive feedback loops mediated by the MDM2 and MDM2 mRNA, respectively, along with the reversible repression of MDM2 caused by Nutlin 3a. For a broader perspective, the reported PBK dynamics of Nutlin 3a are also incorporated. It has been reported that p53 responds to stresses in two ways in terms of concentration to this drug: either it is a sustained (constant) or an oscillatory response. The claimed dosage strategy turned out to be appropriate for sustained p53 response. However, for the induction of oscillations, inhibition of MDM2 is not enough; rather, anti-repression of the p53–MDM2 complex is also needed, which opens new horizons for a new drug design paradigm.

1. Introduction

Cancer is still among the most common causes of death, and it is mostly caused by the stimulation of oncogenes and the inactivity of tumor suppressors. Over the past few years, the tumor suppressor protein p53 has evolved into an important subject in the development of anti-cancer drugs [1,2]. Since the p53 protein was discovered in 1979, scientists have spent a great deal of time and effort investigating its importance in the spread of cancer. Fifty percent of cancer cases have been found to have either p53 protein mutations or inactivation. The protein achieves importance as a result of its involvement in cancer suppression and its capacity to react to a variety of genome-damaging stresses. p53 causes responses such as DNA reparation, apoptosis, and cell death in its wild-type form [3]. When a cell is threatened by a stress (for example, smoking or radioactivity), p53 activates several downstream pathways to guarantee that the cell maintains its normal functional state. Any time the integrity of the genome is in jeopardy, p53 steps in to protect it, earning the nickname “guardian of the genome” [4,5].

The essential function played by p53 in the regulation of many cellular processes necessitates the exact management of its amount at all times. A low steady-state p53 concentration is maintained under normal conditions. p53 is negatively tuned by the murine double minute 2 (MDM2) protein, which also functions as the protein’s E3 ligase [6,7]. An auto-regulatory feedback loop is formed between p53 and MDM2, which regulate each other. MDM2 is the key component in ubiquitination-mediated degradation of p53, whereas when p53 is triggered/activated, it results in the transcription of MDM2 mRNA, which results in the subsequent increase of MDM2 protein [8].

A phosphate ion is bound to the p53 protein by MDM2, causing it to be degraded by the proteasome [9,10,11]. In several tumors, the overexpression of MDM2 is the cause of low p53 levels in many tumors, stopping key cellular processes including DNA reparation, cell cycle stoppage, and apoptosis. As a result, the presence of the protein-to-protein complex between p53 and MDM2 results in the activation and revival of the required strength of wild-type p53, which subsequently revives cell operation in the normal state via p53 moderation. Because of these factors, MDM2 is emerging as a widely utilized therapeutic strategy in cancerous cells [12,13].

The interaction of MDM2 and p53 protein takes place through hydrophobic regions located in specific binding areas [6]. It has been discovered by studying the framework of p53 that certain tiny non-peptide molecules may imitate the complex formation sequence between p53 and MDM2. These substances may interfere with the protein–protein interaction between p53 and MDM2, resulting in an increase in the amount of p53 in the body. Blocking the protein–protein interaction via the use of these molecule inhibitors is becoming a potential treatment approach for human cancers that maintain wild-type p53 [14,15]. In recent years, a large number of small molecule inhibitors have been identified, and several of them have passed clinical and preclinical trials. Nutlin, a class of molecular inhibitors, has a known function of binding to MDM2 such that no genotoxicity can occur. MDM2 binds to Nutlin at the pocket located at the N-terminus, which is the same pocket where p53 attaches [16,17]. Various derivatives of Nutlin such as Nutlin-1, nutlin-2, and nutlin-3 have been discovered by Vassilev et al., and their effects on cancer cells have been discussed in [14]. However, Nutlin-3 is the most effective compound and is commonly used in anti-cancer studies. Nutlin-3a is orally bio-available, and the preclinical data show its good anti-tumor activity, without any toxicity. Nutlin-3a has been shown to restore wild-type p53 functioning, and several additional Nutlin variations have been shown to be helpful in the treatment of tumors harboring defective or mutant p53 [18].

The p53 pathway’s dynamic response is controlled by a complex network of feedback loops. The first investigations were conducted with the goal of monitoring the behaviors at the level of the cell population [18,19]. However, it was later discovered that monitoring the dynamics of a population can obscure the true behavior of individual cells within the population. Therefore, by analyzing a fluorescence-tagged protein, it is possible to deduce the hidden dynamics of a particular cell [20]. The p53–MDM2 loop may show a variety of dynamic response patterns, which vary depending on the input. In general, these patterns may be classified as being oscillatory or sustained [21]. In the event of less severe DNA damage, it has been shown that p53 will start oscillations. Further, the oscillatory behavior may be divided into categories, namely: digital pulses and oscillations with damped or sustained amplitude. It is important to note that the frequency of the pulses is directly related to the amount of DNA damage, but the width of the pulse and its magnitude remain constant [22]. The pulsing p53 gene is often linked with DNA synthesis or apoptosis, respectively. The condition of the DNA is checked after every pulse of about 6 h. In the event that the DNA is fixed, the oscillating p53 is killed, and the cell cycle activity is resumed. The continuous p53 response is triggered as a result of severe DNA damage. The magnitude of the reaction and its width are both in direct proportion to the degree of the damage that has occurred. In this particular instance, the expressed genes result in cell death [22,23]. It is clear that in the event of significant DNA mutation, p53 does not allow the cell to recover in enough time to repair the damage and, instead, kills the cell instantly.

Statistical frameworks enable us to obtain a better understanding of the network structure, develop new ideas, and explore unknown areas. The majority of modeling research on the p53 pathway has focused on the dynamics of p53 and MDM2 [24]. Mathematical models including continuous time, discrete time, and delayed time are used in the literature [25,26,27]. Systems biology (SB) has traditionally been used to analyze and predict the response of complex biological systems. Recently, SB has been utilized as a valuable tool for accurately directing therapeutic treatments in complex biological networks when combined with control theory. While considerable improvements have been achieved in medication delivery to cardiovascular systems [28], blood pressure management [29], tumor chemotherapy [30], anesthetic drug delivery [31], diabetic control [32], Parkinson’s tremor [33], and HIV/AIDS control [34,35], there have also been notable setbacks.

The use of control engineering in cancer therapy is a relatively recent field of study. The primary goal of cancer therapy is to achieve complete remission of malignant cells in the shortest amount of time possible while preserving the patient’s overall health. While chemotherapy, radiation, and surgical procedures are effective treatments, these operations may have a negative impact on the patient’s quality of life [36]. The current research pattern is moving away from traditional methods of analysis and control and toward in silico methods. There is a compelling need to utilize in silico simulations to execute medication design utilizing control theory in order to improve patient outcomes. The control guarantees that the treatment is administered and that the medication is scheduled. Even the most advanced numerical modeling of bio-systems does not always yield an accurate outcome. Consequently, it is necessary to perform the desired task by developing a control system that can cope with the imprecisions and uncertainties inherent in the modeling process. The p53 pathway has been the subject of research in a handful of model-based control methods in the literature. It was shown in [37] how to construct flatness-based control for sustaining a target level of p53 by making use of p53 models and its associated pathways.

It should be noted that none of the control methods described above are fundamentally robust. It is proposed in this article to solve the problem of robustness via the use of the proportional–integral (PI) controller. The block diagram is presented in Figure 1.

2. Mathematical Modeling of p53 Pathway

To describe a dynamic system in the control system paradigm, mathematical modeling is carried out, where a dynamic system can be described in terms of differential equations. Such mathematical models for the p53 pathway have been presented in the literature using different methods including the time-delay-based method and stochastic-based method [21].

Each modeling technique has its own limitations, where the time delay method has to do with real cell proteins and measure the time elapsed during translation and transcription directly. In the stochastic method, rather than taking the instant outcome, the impact of protein levels is quantized, making the system computationally rigorous.

p53 activates MDM2 mRNA by triggering the transcription factor of MDM2. The rate of generation of MDM2 can be controlled by the weighted mean of previous p53 concentration levels. A three-module framework was presented in [38], which would measure the pulse generation of the produced p53, resulting in the breakdown of DNA strands, thus further prompting processes such as DNA reparation, ATM stimulation, and a feedback control loop between p53 and MDM2. Another method of observing the response of p53 on the pathway was presented in [39], which, by generating an input signal, revives p53 and where the subsequent behavior of several upstream processes including ATM, CHK2, and WIPI was observed.

The system dynamics having negative feedback control loops in the p53 pathway was studied in [23], where it was found that the p53–MDM2 complex had sustained oscillations without damping. The research focused on the response of various stresses on the feedback control loop of the p53–MDM2 complex; moreover, the inhibitor was studied, where it was shown to have spiky oscillations. The model inclusively covers all details of the p53 pathway and is suitable enough to be implemented in control system design.

The process of drug development is extremely expensive and time consuming. The drug development process can be enhanced to be more cost effective and accurate. The drug dosage design is also a major problem. For the proposed drug, the dynamics associated with p53 are equally as important as its PK/PD characteristics. Researchers have spent a great deal of effort in physiological-based kinetics (PBK) modeling of numerous drugs. A scheme could be deployed where the PBK characteristics of an anti-tumor drug are integrated with the p53 pathway model. Advances in computational power in recent years have enabled further advances of the in silico models and enhanced the development operation. The complete parametrization of the PBK mathematical model presented by Zhang et al. has been selected for this research.

In the current research, the models presented in the literature by Pszynski and Hunziker are integrated as shown in Figure 2. The p53 pathway dynamics presented in the Hunziker model describe the positive and negative feedback of the mRNA of p53–MDM2. The Pszynski model describes the PBK dynamics of Nutlin. In [37], a new term for the anti-tumor drug Nutlin 3a was added to the Hunziker model and its response investigated.

2.1. p53 Pathway Model by Hunziker

The model concerns four different concentrations (nM): p53 protein, $x_1$; MDM2 mRNA, $x_2$; MDM2, $x_3$; p53–MDM2 complex, $x_4$; these are represented by the following ODEs:

$\begin{array}{l}\frac{d{x}_1}{dt}={\sigma }_{pr}-{\alpha }_i{x}_1-k_{rf}{x}_1{x}_3+k_{rb}{x}_4+{\gamma }_{rd}{x}_4\\ \frac{d{x}_2}{dt}=k_{rt}{x}_1^2-{\beta }_{rd}{x}_2\\ \frac{d{x}_3}{dt}=k_{rl}{x}_2-k_{rf}{x}_1{x}_3+k_{rb}{x}_4+{\alpha }_d{x}_4-K_{B3}{x}_3\\ \frac{d{x}_4}{dt}=k_{rf}{x}_1{x}_3-k_{rb}{x}_4-{\alpha }_d{x}_4-{\gamma }_{rd}{x}_4\end{array}$	(1)
	(2)
	(3)
	(4)

Table 1. Parameters of the p53 model.

Parameters	Description
${\sigma }_{pr}$	Production rate of p53
${\alpha }_i$	MDM2-independent deactivation/degradation of p53
${\alpha }_d$	MDM2-dependent deactivation/degradation of p53
${\gamma }_{rd}$	MDM2 degradation
$k_{rt}$	Transcription of MDM2
${\beta }_{rd}$	Degradation rate of MDM2 mRNA
$k_{rl}$	Translation of MDM2
$k_{rn}$	Nutlin rate constant
$\psi$	Matched disturbance

lists and describes the variables utilized in the p53 model. Here, the symbols (${\alpha }_i$,${\beta }_{rd}$, ${\gamma }_{rd}$, ${\alpha }_d$) depict the degradation rates. The variable ${\alpha }_i$ describes all the actions leading to an MDM2-independent decay of the p53 protein. The variable ${\alpha }_d$ represents the decay of p53 exclusively through MDM2. The MDM2 mRNA decay rate is represented by variable ${\beta }_{rd}$, whereas the variable ${\gamma }_{rd}$ models the decay of the MDM2 protein through the auto-ubiquitination process.

A constant rate is assumed for p53 synthesis and is shown by the subscripted characters, such as the variable ${\sigma }_p$. The transcription rate of MDM2 mRNA is represented by constant $k_{rt}$, whereas the translation rate of the MDM2 protein is represented by $k_{rl}$. The rate constants $k_{rf}$ represent the development of the p53–MDM2 complex, and $k_{rb}$ represents the breakup of the MDM2–p53 complex. Although nearly all variables are constrained, the parameters $k_{rf}$, ${\gamma }_{rd}$, and ${\alpha }_d$ can change due to the environmental parameters of the cell and due to cell–cell fluctuation.

2.2. Modeling PBK Dynamics of Nutlin by Puszynski

Puszynski et al. [25] modeled the response of p53 to the anti-tumor drug Nutlin and determined the extra-cellular concentration and dynamics of the anti-tumor drug with the help of the available pharmacokinetics data of Nutlin oral delivery in mice. It can be noted that due to the complex of Nutlin and plasma proteins, the amount of free Nutlin in plasma is significantly low. The data are incorporated in the following equation:

$$\begin{array}{c}\hfill N_{Qb}=B_p\frac{K_B{N}_{ext}}{1+K_B{N}_{ext}}\end{array}$$

(5)

where $N_{Qb}$ represents the amount of bound Nutlin and $B_p$ and $K_B$ represent the binding site plasma protein and association binding constant having estimated values of $2.86\times {10}^{-4}M$ and $0.085\times {10}^6{M}^{-1}$, respectively. Puszynski et al. took into consideration a uni-compartmental model associated with extra-cellular and intra-cellular parts of a cell, where the total concentration of Nutlin can be calculated as the sum of the concentration of free Nutlin, $N_{ext}$, and concentration of bound Nutlin, $N_{Qb}$.

$$N_{conc}=N_{ext}+N_{Qb}$$

(6)

The extra-cellular free Nutlin $N_{ext}$ can be calculated in terms of $N_{conc}$:

$$\begin{array}{c}\hfill N_{ext}=\frac{-(1+K_B{B}_P-K_B{N}_{conc})+\sqrt{{(1+K_B{B}_P-K_B{N}_{conc})}^2+4K_B{N}_{conc}}}{2K_B}\end{array}$$

(7)

The processes incorporating the intra-cellular free Nutlin can be written as:

$$\begin{array}{c}\hfill \frac{d{n}_c\left(t\right)}{dt}=g_1{N}_{ext}\left(t\right)+k_{diss3}(MD{M}_i\left(t\right)+MD{M}_{pi}\left(t\right)+MD{M}_{pni}\left(t\right))\\ \hfill -k_{B3}{n}_c\left(t\right)(MDM\left(t\right)+MD{M}_p\left(t\right)+MD{M}_{pn}\left(t\right))-e_1{n}_c\left(t\right)\end{array}$$

(8)

where $n_c$ is the concentration of Nutlin in the cell.

Bonding between proteins can also take place in the retina, where it is assumed to be similar to that in plasma. Therefore, the four assumptions are as follows:

• Drug dissemination takes place in a single compartment.
• The elimination of only the free Nutlin takes place.
• The elimination is linear.
• The binding proteins are in a quasi-steady state.

This results in the following pharmacokinetics equation:

$$\begin{array}{c}\hfill \frac{d{N}_{conc}}{dt}=p_m{D}_r\left(t\right)-{\alpha }_2{N}_{ext}\left(N_{conc}\right),N_{conc}\left(t_0\right)=0\end{array}$$

(9)

where $p_m$ is the conversion factor from $mgK{g}^{-1}$ to $moles$, $D_r$ is the rate of drug dosage, $t_0$ is the initial time of delivery, and $N_{conc}$ is according to Equation (6). When dealing with oral delivery, the gastro-enteric release is assumed to be exponential; hence, for a single dose of operation at $t=t_0$, the equation can be re-written as the following equation:

$$\begin{array}{c}\hfill \frac{d{N}_{conc}}{dt}=p_m{D}_{dose}{\alpha }_1{e}^{-{\alpha }_1(t-t_0)}-{\alpha }_2{N}_{ext}\left(N_{conc}\right),N_{conc}\left(t_0\right)=0\end{array}$$

(10)

where $D_{dose}$ is the dose (in $mgK{g}^{-1}$).

The mathematical model presented in this manuscript considers the negative feedback loop involving p53 and its inhibitor, Mdm2, at the core of the p53 pathway, and uses it to examine the effect of Nutlin3-a to enhance the activity of the p53 protein. Therefore, this study is applicable to those cancer types where the cancer is caused by the inactivation of the p53 protein due to the overexpression of MDM2. The other cancer types that are caused by the mutation of the p53 gene or involve any other signal transduction networks will have completely different mathematical models involving proto-oncogenes as state variables and are outside of the scope of this study.

3. Control Methodology

A nested feedback control technique comprising two loops, as shown in Figure 3, is implemented to produce a sustained response of the p53 pathway. The outer loop consists of the p53–MDM2 pathway and its controller. The purpose of this control loop is to determine the required amount of Nutlin ($n_{ref}$) in the pathway to achieve the desired p53 response. This provides us the reference dosage of Nutlin that must be present inside the cell to achieve the desired p53 response. However, as discussed previously in Equation (6), the total extracellular Nutlin concentration is decomposed into two portions. One major portion is the concentration of blood-plasma-bound Nutlin, $N_{Qb}$, and the other is the extra-cellular free Nutlin, $N_{ext}$. The later one is the only available portion of Nutlin that can be transferred inside the cell.

For the inner control loop, the equation becomes:

$$\begin{array}{c}\hfill u_{nd}\left(t\right)=P+I+D\end{array}$$

(11)

$$P=k_p{e}_{nc}\left(t\right)$$

(12)

$$I=k_i{\int }_0^t{e}_{nc}\left(\tau \right)d\tau$$

(13)

$$D=k_d\frac{d{e}_{nc}\left(t\right))}{dt}$$

(14)

$$u_{nd}\left(t\right)=k_p{e}_{nc}\left(t\right)+k_i{\int }_0^t{e}_{nc}\left(\tau \right)d\tau +k_d\frac{d{e}_{nc}\left(t\right)}{dt}$$

(15)

The outer loop equation becomes:

$$n_{ref}\left(t\right)=k_{p}e\left(t\right)+k_i{\int }_0^{t}e\left(\tau \right)d\tau +k_d\frac{de\left(t\right)}{dt}$$

(16)

Therefore, to maintain the reference Nutlin dosage $n_{ref}$ inside the cell, another feedback loop was introduced. The inner feedback control loop is designed to consider the Nutlin dynamics of the PBK model. In this cascaded control arrangement, both the inner and outer loops are running simultaneously, where the outer loop generates the reference Nutlin dosage and the inner loop ensures that this dosage is received inside the cell considering the cellular dynamics of the drug. It is pertinent to mention that in this cascade arrangement, the inner loop should be reasonably fast as compared to the outer loop, so that the inner loop can effectively track all the changes in the reference. We require a simple and fast controller with a relatively easy implementation, especially in vivo. Therefore, we implemented the proportional–integral–derivative (PID) controller in both feedback loops. The outer PID controller is the function of the error between the desired state vector $X_{ref}$ and the output state vector $X\left(t\right)$. In the same fashion, the inner feedback loop PID controller is a function of the error between $n_{ref}$ and the Nutlin present in the cell $n_c$.

4. Results and Discussions

The results are based on the PID approach in the domain of control systems. The drug concentration is calculated through a two-fold system of negative feedback loops, as shown in Figure 3. The first loop is targeted towards the analysis of MDM2 concentrations, which are directly involved in the computations of the reference Nutlin 3a dose. The reference dose is obtained through the outer loop PID controller. In order to maintain appropriate dose concentrations, the second loop comes into play to provide the information about the Nutlin 3a levels already present within the body. With the inclusion of this final step, the dose concentration becomes the function of the error between the reference dosage $n_{ref}$ and the actual Nutlin concentration present in the cells.

As discussed earlier, there are two ways in which p53 can be triggered depending on the level of stress. If the stress levels are not very significant, then the response of the p53 pathway is oscillatory, resulting in several subsequent downstream reactions that stop the call cycle and mend the cell. If the damage is repaired, then the p53 oscillations are halted and the cell cycle continues in its normal fashion. However, when the damage to the cell is quite significant, the p53 pathway results in sustained oscillations, resulting in apoptosis.

However, it has been reported that p53 is not present or mutated in 50% of cancer cases. This is due to the over-expression of MDM2, which binds with p53 through a negative feedback, resulting in the degradation of p53 through the ubiquitination process.

The problem could be countered in such way that another molecule can be introduced that binds with MDM2 in the same way as p53. Nutlin, a drug molecule, possesses similar characteristics. The process of achieving the required response of the p53 pathway due to Nutlin can be divided into two strategies. The first one is sustained response, while the second one is incorporating the p53 response.

In order to repair the cells, when the damage to the DNA is miniscule, p53 should respond in an oscillatory manner. The degeneration rate of MDM2 and p53–MDM2 increases and decreases, respectively, in response to DNA damage.

The effectiveness of the suggested control method was assessed using closed loop simulations, as shown in Figure 4, Figure 5, Figure 6 and Figure 7. The concentrations of all state variables are clearly approaching their desired equilibrium values. The convergence time was less than an hour, which may be too short given the substantial overshoot of the Nutlin dose. This can be decreased by modifying the gains of the controllers according to therapeutic needs.

5. Conclusions

In order to suppress the development of the p53–MDM2 structure, a model that is founded on proportional–integral–derivative (PID) control was used in conjunction with the PBK model of Nutlin 3a, which led to the calculation of its dose concentration based on the damage of p53. A two-step process was devised, which first calculates the reference Nutlin 3a level that should be available in the body to counter the MDM2 levels, while the second step comprises feedback from the body itself to enable the calculation of the differential amount of the required Nutlin 3a. Based on in silico trials, sustained p53 response was validated. For this study, the MDM2 degradation constant was analyzed, and it was reported that the anti-repression of the p53–MDM2 complex is also required. By and large, it was shown that the apoptosis of the affected cells is possible through the injection of the right amount of Nutlin 3a, which also conforms to the safe levels of the substance within the body.

However, repairing of the damaged cells is not possible as it requires breaking the p53–MDM2 bond. To counter this, new drug structures need to be explored. Moreover, future work could be extended in the direction of developing an enhanced mathematical model that keeps in view the complex networks of proto-oncogenes having different signaling pathways, along with including the possibility of Nutlin administration in the early stages of cancer, when p53 is still functional.