Mathematical modelling of Her2 (ErbB2) PI3K/AKT signalling pathways during breast carcinogenesis to include PTPD2

1.   Introduction

Breast cancer is one of the most common cancers, a main cause for mortality due to cancer in women worldwide ^[1]. Central to many forms of breast cancer is the altered activities of the ErbB family members which are receptor tyrosine kinases, all structurally related to the epidermal growth factor receptor (EGFR) including Her1 (EGFR, ErbB1), Her2 (Neu, ErbB2), Her3 (ErbB3), and Her4 (ErbB4) ^[2]. About 25% of breast cancer patients display an overexpression of ErbB2, which plays a pathogenic role in breast cancer ^[3]. ErbB2 overexpression at the early stages of the disease is linked to a poor prognosis and is related to the development of invasive tumors ^[4]. Therefore, understanding of the regulation of ErbB2 signalling is essential for the treatment of breast cancer.

ErbB2 (HER2/Neu) is the only receptor in the EGFR family that has no specific ligand ^[5]. Where ErbB2 is overexpressed, it will self-activate or can form a heterodimer with ErbB1, ErbB3 or ErbB4. Following such events, ErbB2 then undergoes autophosphorylation, producing a phosphotyrosine site that recruits and activates signalling complexes. This leads to the activation of the effector pathway, mainly the RAS/mitogen activated kinase (MAPK) and the phosphoinositide 3-kinase (PI3K)/protein kinase B (AKT) signalling pathways ^[6]. Mathematical modelling is a popular method to study signal pathways, especially for complex signal networks ^[7,8]. To better understand how the many ErbB1-linked signalling components function in different contexts, mathematical modelling is used to provide new insights. Thus, Schoeberl et al. established a computational model for EGF receptor-induced MAPK kinase pathway, suggesting that the signal efficacy is primarily dependent on the initial velocity of receptor activation ^[9]. Others ^[10] have generated a computer simulation model of the HRG (Heregulin Beta1) induced ErbB4 signalling and model simulations indicate that the HRG pathway is regulated by PP2A (protein phosphatase 2A) and the interaction between RAF and AKT.

The same PI3K / AKT and MAPK pathways and crosstalk mechanisms were modeled using the parameter estimation decomposition method by mixed-function Petri nets ^[11], and supported a new hypothesized interacting mechanism between the Akt pathway to the MAPK pathway. A different approach based on ODEs (ordinary differential equations) modelling was used for the HER2/3-AKT signalling pathway in order to study cancers that overexpress HER2 ^[12]. This model helped to understand evaluate the efficacies of several therapies and further proposed a potential “steering-based” treatment scheme for cancers that overexpress HER2. The importance of signalling pathway modelling was illustrated in a ErbB2/3 network model which was used to design a treatment regimen for ErbB2-positive cancer by using quantitative logic. This indicated that a combination of inhibitors which targeted ErbB3 and ErbB2 is more effective than targeting AKT and the mitogen-activated protein kinase MEK ^[13]. Recently, the PI3K-Akt and NF-κB signaling pathways have been found to play a key role in the progression of prostate cancer. Shankar et al. used a complex systems biology approach that focused on a multi-layered, hierarchical paradigm to model these pathways modularized to identify the components that played a coordinating role in the two pathways ^[14].

These studies have modeled several aspects of the ErbB2 signalling network but on-going research demands that such models need to be constantly updated, especially if the clinical treatment of breast cancer has been optimized. Clinical treatment of ErbB2-positive breast cancer is mainly based on the drugs that specifically target the ErbB2 protein, such as trastuzumab (herceptin). Inhibiting ErbB2 activity results in tumor cells showing a degree of apoptosis but surviving cells develop drug resistance in about 80% of patients ^[15]. A new positive regulator of ErbB2 signalling network, PTPD2 has been found to act via PI3K / AKT signalling ^[16] which crucially has been associated with resistance to drug therapies ^[17]. Given this key new finding, we here derive a mathematical model for the PI3K / AKT signalling pathway which incorporates PTPD2 axis to explore the treatment of ErbB2-positive breast cancer and resistance to its treatment.

2.   Materials and method

2.1.   Basic structure of the model

In this model, it is assumed that the activation of the PI3K / AKT signalling pathway and PTPD2 arises following phosphorylation of ErbB2. This is shown in details in Figure 1 which also highlights the lack of ligand binding. The phospholipase D2 (PLD2) /PTPD2 pathway is described in the upper right of Figure 1. PLD2 is an enzyme involved in the processing of phosphatidylcholine (PC) to produce phosphatidic acid (PA) and choline (C) ^[18]. Following activation, PLD2 is recruited to the plasma membrane where it binds to PC which is hydrolyzed to PA ^[19]. PA binding to PTPD2 enhances its phosphatase activity ^[16] to influence in the downstream pathway of ErbB2.

Figure 1.  A schematic description of the ErbB2 signalling pathway.

DownLoad: Full-Size Img PowerPoint

The PI3K/AKT pathway is described by the down left side of Figure 1. Orthodox ErbB2 pathways include RAS/MAPK and PI3K/AKT ^[20]. However, Ramesh et al. found that AKT signalling inhibited PTPD2, whilst MAP and ERK, did not change significantly ^[16]. Therefore, we here concentrate only on modelling the PI3K/AKT signalling pathway, without the RAS/MARK component to make the model more concise. Model of PI3K/AKT pathway involves the following consecutive steps: 1) The activated phosphoinositide3-kinase (PI3K) dimerizes with phosphatidylinositol-4, 5-biphosphate (PIP2) to form a dimer PIP2-PI3K. 2) PIP2 is autophosphorylated on PIP2-PI3K that decomposes into PI3K and phosphatidylinositol-3, 4, 5-triphosphate (PIP3) ^[21]. 3) Phosphatase and tensin homolog (PTEN) removes the phosphorus from PIP3, by binding with PIP3 and forming a complex PIP3-PTEN, and thus changes PIP3 into PIP2. 4) PIP3 and AKT combine into a complex PIP3-AKT, on which AKT can phosphorylate to ppAKT ^[22]. 5) PP2A is an enzyme that dephosphorylates ppAKT to AKT by binding to ppAKT to form PP2A-ppAKT, on which ppAKT is dephosphorylated to produce PP2A and AKT.

How AKT and PTPD2 interact has not been documented, but can be inferred from the experimental results of Ramesh et al. that phosphorylated AKT decreases the inhibition of PTPD2 group without a change in phosphorylated ERK ^[16]. The PI3K/AKT pathway interacts with the RAF/MEK/ERK pathway: PIP3 can activate MEK through PDK1, which affects ERK phosphorylation ^[11]. PTP SHP-1 which belongs to Classical PTPs with PTPD2 can be combined with both PA and PIP3 ^[23]. We suggest that PTPD2-PA promotes the phosphorylation of AKT by binding to PIP3-AKT to form a PIP3-AKT-PTPD2-PA complex. AKT phosphorylation of this complex decomposes into ppAKT, PTPD2, PIP3 and PA, which has a higher activation rate than that of PIP3-AKT.

2.2.   Model equations

The model uses 22 ODEs to simulate the biochemical reaction of the PI3K/AKT signalling pathway and the PLD2/PTPD2 pathway, and the relation between them after activation (Figure 1). The first seven formulas (Eqs (1)–(8)) describe the changes in the concentrations of each component of the PLD2/PTPD2 signalling pathway with time. The remaining 14 formulas (Eqs (9)–(22)) describe the variation in the concentrations of each component of the PTPD2-affected PI3K/AKT signalling pathway. The hydrolysis process of phosphatidylcholine (PC) is based on interfacial kinetics and interfacial binding is considered ^[24]. According to these theoretical considerations the following equations are presented:

Based on Adi et al.’s mathematical model of phosphorylation AKT ^[25] and the hypothesis that PTPD2 promotes AKT phosphorylation by forming a PIP3-AKT-PTPD2-PA complex, our model utilizes the following 14 ODEs to simulate the PI3K/AKT signalling pathway affected by PTPD2.

● PLD2(t), PL(t), PC(t), PA(t), C(t) and PTPD2(t) represent concentrations of proteins of PLD2, PL, PC, PA and PTPD2, respectively.

● PL~PLD2(t), PL~PLD2~PC(t) and PA~PTPD2 represent concentrations of these bi-molecular or tri-molecular complexes, respectively.

● PI3T(t), PIP2(t), PIP3(t), PTEN(t), AKT(t), ppAKT(t) and PP2A(t) represent concentrations of proteins of PI3K, PIP2, PIP3, PTEN, ppAKT, AKT and PP2A respectively.

● PIP2~PI3K(t), PIP3~PTEN(t), PIP3(t)~AKT(t), PIP3~AKT~PTPD2~PA(t), PIP3~PTPD2~PA(t) and PP2A~ppAKT(t) represent concentrations of these bi-molecular tri-molecular or four-molecular complexes respectively.

Equations (1)–(22) describe the variation in concentrations of proteins and complexes involved in the PLD2/PTPD2 and PI3K/AKT pathways with time, respectively. For example, $\frac{d}{dt}PLD2\left(t\right)$ is the variation of PLD2 with time as demonstrated in Eq (1), ${b}_{1}$ represents the binding rate constant of PL and PLD2, ${d}_{1}$ represents the dissociation rate constant of complex PL~PLD2. The definitions of all parameters in Eqs (1)–(22) are described in details in Table 1.

Table 1.  Definitions of model parameters used in the model.

Parameter	Description
${\mathit{\boldsymbol{b}}}_{1}$	The binding rate constant of PL and PLD2
${\mathit{\boldsymbol{d}}}_{1}$	The dissociation rate constant of complex PL~PLD2
${\mathit{\boldsymbol{b}}}_{2}$	The binding rate constant of PL~PLD2 and PC
${\mathit{\boldsymbol{d}}}_{2}$	The dissociation rate constant of complex PL~PLD2~PC
${\mathit{\boldsymbol{k}}}_{1}$	The hydrolysis rate constant of PC
${\mathit{\boldsymbol{b}}}_{3}$	The binding rate constant of PTPD2 and PA
${\mathit{\boldsymbol{d}}}_{3}$	The dissociation rate constant of complex PA~PTPD2
${\mathit{\boldsymbol{b}}}_{4}$	The binding rate constant of PIP2 and PI3K
${\mathit{\boldsymbol{d}}}_{4}$	The dissociation rate constant of complex PIP2~PI3K
${\mathit{\boldsymbol{k}}}_{2}$	The phosphorylation rate constant of PIP2 on complex PIP2~PI3K
${\mathit{\boldsymbol{b}}}_{5}$	The binding rate constant of PIP3 and PTEN
${\mathit{\boldsymbol{d}}}_{5}$	The dissociation rate constant of complex PIP3~PTEN
${\mathit{\boldsymbol{k}}}_{3}$	The dephosphorylation rate constant of PIP3 on complex PIP3~PTEN
${\mathit{\boldsymbol{b}}}_{6}$	The binding rate constant of PIP3 and AKT
${\mathit{\boldsymbol{d}}}_{6}$	The dissociation rate constant of complex PIP3~AKT
${\mathit{\boldsymbol{b}}}_{7}$	The binding rate constant of PIP3~AKT and PA~PTPD2
${\mathit{\boldsymbol{d}}}_{7}$	The dissociation rate constant of complex PIP3~AKT~PTPD2~PA
${\mathit{\boldsymbol{k}}}_{4}$	The phosphorylation rate constant of AKT on complex PIP3~AKT~PTPD2~PA
${\mathit{\boldsymbol{k}}}_{5}$	The dissociation rate constant of complex PIP3~PTPD2~PA
${\mathit{\boldsymbol{b}}}_{8}$	The binding rate constant of PP2A and ppAKT
${\mathit{\boldsymbol{d}}}_{8}$	The dissociation rate constant of complex PP2A~ppAKT
${\mathit{\boldsymbol{k}}}_{6}$	The dephosphorylation rate constant of ppAKT on complex PP2A~ppAKT

---

 | Show Table

DownLoad: CSV

3.   Results

The values of parameters used in the model equations are described in Table 2. Following the previous kinetic analysis of PLD by Henage et al. ^[26], the equation parameter values of PC axis are obtained through conversion and calculation. The equation parameter values of the PI3K / AKT signalling pathway are derived from the computational model of the AKT pathway developped by Wan et al. ^[27]. In addition, the parameter values of the equations related to PTPD2 are obtained by fitting the experimental data.

Table 2.  Values of parameters used in the model equations.

Parameters	Values	Units	Comments
${\mathit{\boldsymbol{b}}}_{1}$	1.0e-6	$\mathrm{n}{M}^{-1}{s}^{-1}$	[26]
${\mathit{\boldsymbol{d}}}_{1}$	1.0+0	${s}^{-1}$	[26]
${\mathit{\boldsymbol{b}}}_{2}$	1.0e-3	$\mathrm{n}{M}^{-1}{s}^{-1}$	[27]
${\mathit{\boldsymbol{d}}}_{2}$	1.0e-1	${s}^{-1}$	[27]
${\mathit{\boldsymbol{k}}}_{1}$	3.6e+3	${s}^{-1}$	[28]
${\mathit{\boldsymbol{b}}}_{3}$	5.0e-4	$\mathrm{n}{M}^{-1}{s}^{-1}$	Fitting
${\mathit{\boldsymbol{d}}}_{3}$	1.0e-1	${s}^{-1}$	Fitting
${\mathit{\boldsymbol{b}}}_{4}$	5.0e-6	$\mathrm{n}{M}^{-1}{s}^{-1}$	[29]
${\mathit{\boldsymbol{d}}}_{4}$	1.0e-1	${s}^{-1}$	[29]
${\mathit{\boldsymbol{k}}}_{2}$	2.0e-1	${s}^{-1}$	[29]
${\mathit{\boldsymbol{b}}}_{5}$	5.0e-6	$\mathrm{n}{M}^{-1}{s}^{-1}$	[29]
${\mathit{\boldsymbol{d}}}_{5}$	1.0e-1	${s}^{-1}$	[29]
${\mathit{\boldsymbol{k}}}_{3}$	1.0e-1	${s}^{-1}$	[29]
${\mathit{\boldsymbol{b}}}_{6}$	2.6e-4	$\mathrm{n}{M}^{-1}{s}^{-1}$	[29]
${\mathit{\boldsymbol{d}}}_{6}$	1.0e-1	${s}^{-1}$	[29]
${\mathit{\boldsymbol{b}}}_{7}$	5.0e-6	$\mathrm{n}{M}^{-1}{s}^{-1}$	Fitting
${\mathit{\boldsymbol{d}}}_{7}$	1.0e-1	${s}^{-1}$	Fitting
${\mathit{\boldsymbol{k}}}_{4}$	3.0e+0	${s}^{-1}$	Assumption
${\mathit{\boldsymbol{k}}}_{5}$	2.0e-1	${s}^{-1}$	Assumption
${\mathit{\boldsymbol{b}}}_{8}$	1.7e-6	$\mathrm{n}{M}^{-1}{s}^{-1}$	[29]
${\mathit{\boldsymbol{d}}}_{8}$	1.0e-1	${s}^{-1}$	[29]
${\mathit{\boldsymbol{k}}}_{6}$	1.5e+0	${s}^{-1}$	[29]

---

 | Show Table

DownLoad: CSV

The initial values of proteins in the model are listed in Table 3. They are amended to match the simulation results with the experimental data. In order to obtain the change of components concentration over time observed in the experimental data, Image J was used to process the protein electrophoretogram of time course experiment of the ERBB2 signaling pathway medium activation in the work of Ramesh et al. ^[16]. Then the relative concentrations of PTPD2 and ppAKT before and after PTPD2 knockout were obtained respectively. These data were used to serve as initial levels of corresponding components based on which we ran our simulations. The model equations were simplified in the simulation by assuming a quasi-steady-state where the concentration of enzyme - base complex remained constant. The simulation was performed using the Matlab software package (R2015b, Mathworks, Natick, USA), and the mathematical equations of the model were solved by the ode45 solver.

Table 3.  Initial concentrations of proteins used in the simulation.

Protein	Value(𝐧nM)	Protein	Value(nM)
PL	2.0e+5	PIP2	1.0e+4
PLD2	1.0e+2	PIP3	0
PC	4.0e+2	PTEN	3.5e+3
PA	1.0e+2	AKT	7e+3
C	1.0e+1	ppAKT	0
PTPD2	2.0e+3/2.4e+2	PP2A	4e+3
PI3K	1.0e+2

---

 | Show Table

DownLoad: CSV

Figure 2 shows the temporal variation in the normalized concentrations of ppAKT before and after PTPD2 knockdown, respectively, based on model simulation results, together with changes of normalized corresponding experimental data. The initial concentration of PTPD2 is set as 2000 nM and 240 nM before and after PTPD2 knockdown respectively. Figure 3 and Figure 4 demonstrate the variation of ppAKT concentration with different inintal concentration of PTPD2 and PLD2, respectively, aiming to investigate how the variations of PTPD2 and PLD2 affect ppAKT. The initial concentration of PTPD2 is set as 2000nM in Figure 4. Figure 5 presents the temporal variations in concentration of related components in the PI3K/AKT signalling pathway including PTPD2. The initial concentration of PTPD2 is set as 2000 nM. In order to observe the influence of PI3K and PTEN on ppAKT, Figure 6 and Figure 7 illustrate the temporal variations in ppAKT concentration with different initial concentrations of PI3K and PTEN, respectively.

Figure 2.  Model simulation of changes in normalized concentrations of ppAKT in the PI3K/AKT signalling pathway containing PTPD2, together with variations normalized corresponding experimental data before and after PTPD2 knockdown (The initial concentration of PTPD2 was set as 2000nM and 240nM, respectively).

DownLoad: Full-Size Img PowerPoint

Figure 3.  Model simulation of changes in ppAKT concentration with different initial concentrations of PTPD2.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Model simulation of changes in ppAKT concentration with different initial concentrations of PLD2 (The initial concentration of PTPD2 was set as 2000 nM).

DownLoad: Full-Size Img PowerPoint

Figure 5.  Model simulation of changes in components’ concentration in the PI3K/AKT signalling pathway including PTPD2 (The initial concentration of PTPD2 was set as 2000 nM).

DownLoad: Full-Size Img PowerPoint

Figure 6.  Model simulation of changes of ppAKT concentration with different initial concentrations of PI3K (The initial concentration of PTPD2 was set as 2000 nM).

DownLoad: Full-Size Img PowerPoint

Figure 7.  Model simulation of changes in ppAKT concentration with different initial concentrations of PTEN (The initial concentration of PTPD2 was set as 2000 nM).

DownLoad: Full-Size Img PowerPoint

4.   Discussion

By modelling the PI3K/AKT signalling pathway and the newly discovered PTPD2 axis, the variation in the concentration of each component in the new PI3K/AKT signalling pathway over time after ErbB2 activation can be obtained. Further, the effect of other phosphatases on AKT phosphorylation can be investigated as well.

These following results substantiate the validity of the established mathematical model. On the one hand, experimental results of Ramesh et al. ^[16] observations indicated that AKT was rapidly phosphorylated to ppAKT (within 5 minutes) under the acute activation of ErbB2 and that the inhibition of PTPD2 resulted in reduced AKT phosphorylation. Similar to the experiment by Ramesh et al. ^[16], after reducing the initial concentration of PTPD2 by about one-fifth, it was observed that the steady-state concentration of ppAKT also decreased by about one third. From Figure 2, we can see that the difference between the model simulation results and the experimental data is relatively small. As demonstrated in Figure 3 and Figure 4, the decrease in the initial concentration of PTPD2 and PLD2 delays the time to reach the peak for ppAKT or reduces the steady-state value of ppAKT, thereby reducing the phosphorylation level of AKT. This concentration effect is important as high levels of ppAKT drive the appearance of the multiacinar phenotype which has consequences in terms of carcinogenesis and treatment ^[30,31]. These results match the experimental conclusion of Ramesh et al. ^[16] that inhibition of PTPD2 and PLD2 can reduce the multiacinar phenotype of mammary epithelial cells. On the other hand, as shown in Figure 5, the simulation results show that the activaton of the PI3K/AKT pathway leads to an initial decreases and increases in PIP2 and AKT (PIP3 and ppAKT), respectively, and then a steady state after about 400 seconds. These results are consistent with the accepted modelling results by Adi et al. ^[25] of PI3K / AKT pathway without the new regulatory factors.

Furthermore, high concentrations of ppAKT have been proven to lead to enhanced anti-apoptotic properties of cells. Therefore, the model was used to simulate the changes of ppAKT concentration over time at different initial concentrations of PTPD2 and PLD2. The simulation results show that the lower the initial concentration of PTPD2 and PLD2 was, the lower the time-integrated concentration of ppAKT (Figures 3 and 4). Then, the effect of the initial concentrations of other enzymes in the signalling pathway on ppAKT concentration was also studied. It is found that drop in the initial concentration of PI3K leads to decrease in the concentration of ppAKT. Whereas, the decrease in the initial concentration of PTEN causes the increase in the concentration of ppAKT (Figures 6 and 7). Therefore, high concentrations of PI3K, PTPD2 or PLD2 can rapidly drive ppAKT accumulation. The steady levels of ppAKT that are achieved appear to be influenced by PI3K and PTPD2 concentrations.

5.   Conclusions

ErbB2 overexpression plays an important pathogenic role in breast cancer. In ErbB2 positive breast cancer, the ErbB2 signalling network affects cell proliferation growth and survival ^[32]. Analysis of this signalling network by mathematical modelling will therefore be useful in understanding how components in this pathway contribute towards carcinogenic events. The PI3K/AKT signalling pathways have been previously mathematically modeled but fail to include the newly discovered the role of PTPD2 component. In this context, we derived a new mathematical model incorporating PTPD2 in the PI3K/AKT signalling pathway to explore the couplling mechanism between two pathways.

In our model, we assumed that PTPD2 acts on the AKT signalling pathway by combining PA with PTPD2 and participates in AKT phosphorylation with PIP3. Based on this, we established a mathematical model including both PI3K/AKT and PLD2/PTPD2 pathways, and conducted simulation based on this model. Crucially, the simulation results were consistent with the experimental results of Ramesh et al. ^[16]. The model also partly explains why inhibiting PTPD2 can attenuate the multiacinar phenotype of mammary epithelial cells. We also studied the effect of the initial concentration of phosphatase on the change in ppAKT concentration. Taken together, these observations suggest therapeutic approaches for erbb2-positive breast cancer that is resistant to ErbB2 targeted therapy based on inhibitors for PI3K, PTPD2 or PLD2.

Our model does not include the downstream pathway of ppAKT, nor does it establish a connection between the pathways or effects at the cellular level. Here, we have only demonstrated that the over-activation of PI3K and PTPD2 induces over-activation of ppAKT based on model simulations. In the future, future work is preferred to establish a more complete mathematical model of ErbB2 signalling network, exploring the relationship between the signalling pathway and breast cells.

Acknowledgments

This work was supported by National Key R & D Program of China [grant 2018YFB1305400]; the National Natural Science Foundation of China [grant numbers 61673246, 81301294, U1806202 and 61533011]; and the Research and Development Program of Jinan [grant number 201907064].

Conflict of interest

We declare that we have no competing interests.