Mathematical modelling of the role of GADD45β in the pathogenesis of multiple myeloma — TEI-tagged view

Yao ZhangpersName, Changqing ZhenpersName, Qing YangpersName, Bing JipersName (2020date). “Mathematical modelling of the role of GADD45β in the pathogenesis of multiple myeloma.” Royal Society Open Sciencetitle·j.

Object of study — the interaction of NF-κBterm·dom and JNKterm·dom signalling through GADD45βterm·dom and MKK7term·dom in multiple myelomaterm·dom cells and bone

Research question — how GADD45βterm·dom, downstream of NF-κBterm·dom, promotes MMterm·dom cell survival through the JNKterm·dom pathway, a mechanism earlier models ignored

Contribution — an ODEterm·mth model coupling bone cells, MMterm·dom cells and the NF-κBterm·dom and JNKterm·dom pathways, used to compare GADD45β- and NF-κB-targeted therapies

Full text — tagged

Mathematical modelling of the role of GADD45βterm·dom in the pathogenesis of multiple myelomaterm·dom Abstract Multiple myelomaterm·dom (MMterm·dom) is an incurable disease with relatively high morbidity and mortality rates. Great efforts were made to develop nuclear factor-kappa B (NF-κBterm·dom)-targeted therapies against MMterm·dom disease. However, these treatments influence MMterm·dom cells as well as normal cells, inevitably causing serious side effects. Further research showed that NF-κBterm·dom signalling promotes the survival of MMterm·dom cells by interacting with JNKterm·dom signalling through growth arrest and DNA damage-inducible beta (GADD45βterm·dom), the downstream module of NF-κBterm·dom signalling. The GADD45β-targeted intervention was suggested to be an effective and MMterm·dom cell-specific treatment. However, the underlying mechanism through which GADD45βterm·dom promotes the survival of MMterm·dom cells is usually ignored in the previous models. A mathematical model of MMterm·dom is built in this paper to investigate how NF-κBterm·dom signalling acts along with JNKterm·dom signalling through GADD45βterm·dom and MKK7term·dom to promote the survival of MMterm·dom cells. The model cannot only mimic the variations in bone cells, the bone volumeterm·dom and MMterm·dom cells with time, but it can also examine how the NF-κBterm·dom pathway acts with the JNKterm·dom pathway to promote the development of MMterm·dom cells. In addition, the model also investigates the efficacies of GADD45βterm·dom - and NF-κB-targeted treatments, suggesting that GADD45β-targeted therapy is more effective but has no apparent side effects. The simulation results match the experimental observations. It is anticipated that this model could be employed as a useful tool to initially investigate and even explore potential therapies involving the NF-κBterm·dom and JNKterm·dom pathways in the future.

Keywords: mathematical model, multiple myelomaterm·dom, GADD45βterm·dom|MKK7term·dom complex, the NF-κBterm·dom pathway Introduction Multiple myelomaterm·dom (MMterm·dom) is a neoplastic plasma cell disease with high morbidity and mortality rates. MMterm·dom reportedly accounts for approximately 10–15% of haematologic malignancies in the world [ 1 , 2 ]. Almost 6–7 individuals out of every 100 000 people in the world are diagnosed with MMterm·dom each year [ 2 ]. In 2015, 28 850 patients were diagnosed with MMterm·dom, and the number of deaths due to this illness was more than 11 000 in the United States [ 3 ]. However, MMterm·dom is still incurable and there is a great need to develop more effective therapeutic strategies.

Nuclear factor-kappa B (NF-κBterm·dom) was reported to be abnormally activated in many types of cancer cells, including MMterm·dom in the 1990s [ 4 ], and it inhibits the apoptosisterm·key of cancer cells by upregulating anti-apoptotic genes [ 4 , 5 ]. A series of studies have been performed to investigate the underlying mechanism of NF-κBterm·dom in the pathogenesis of MMterm·dom, with the aim of exploring NF-κB-targeted therapies [ 6 – 11 ]. However, to date, no appropriate NF-κB-targeted treatments have been established due to their serious side effects [ 4 ]. The inhibitors of NF-κBterm·dom similarly affect the NF-κBterm·dom pathway in MMterm·dom and normal cells leading to serious damage. This reflects the important roles of NF-κBterm·dom in cell survival, inflammatory and immune responses and other issues [ 12 ].

Growth arrest and DNA damage-inducible beta (GADD45βterm·dom) forms part of the downstream module of the NF-κBterm·dom pathway and is essential to the survival of MMterm·dom cells [ 4 , 13 ]. GADD45βterm·dom binds to and inhibits MAP kinase kinase 7 (MKK7term·dom) to reduce the phosphorylation of c-Jun N-terminal kinase (JNKterm·dom) pathway which usually acts to ensure the normal apoptosisterm·key of cell [ 14 ]. Therefore, the overexpression of GADD45βterm·dom due to an abnormally activated NF-κBterm·dom pathway in MMterm·dom cells can inhibit the apoptosisterm·key of MMterm·dom cells by suppressing the JNKterm·dom pathway [ 4 , 14 , 15 ]. Hence, GADD45βterm·dom|MKK7-targeted therapeutic strategies were suggested as a potential way to kill MMterm·dom cells effectively, and crucially, these may not have side effects on normal cells [ 4 , 16 ].

A series of mathematical models have been constructed to study MMterm·dom based on experimental findings [ 17 – 23 ]. These models have demonstrated great potential to improve our understanding of the complicated pathogenesis of the disease. However, the underlying mechanism through which GADD45βterm·dom promotes the survival of MMterm·dom cells as downstream modules of the NF-κBterm·dom pathway was not included. Therefore, a mathematical model of MMterm·dom is built in this paper by considering the NF-κBterm·dom and JNKterm·dom pathways and their interaction through GADD45βterm·dom and MKK7term·dom. The model cannot only mimic the variations in bone cells, the bone volumeterm·dom and MMterm·dom cells with time, but it can also examine how the NF-κBterm·dom pathway acts with the JNKterm·dom pathway to promote the development of MMterm·dom cells.

Model development Basic structure of the model Figure 1 describes the schematic diagram of the constructed mathematical model in the paper. It demonstrates the mechanisms of the NF-κBterm·dom and JNKterm·dom pathways and shows how the interaction between these two pathways promotes the survival of MMterm·dom cells. It should be noted that the coupling between MMterm·dom cells and bone cells (osteoclastic and osteoblastic lineages) is not included in figure 1 , but it is included in the electronic supplementary material, Figure A1 in Appendix A for convenience, since these cells have been previously studied (e.g. [ 23 ]). Figure 1 is made up of two parts: part A and part B, which involve the NF-κBterm·dom and JNKterm·dom pathways, respectively.

The schematic description of the mechanisms of the NF-κBterm·dom and JNKterm·dom pathways, and their interaction.

Part A shows how the NF-κBterm·dom pathway is activated and how it then upregulates anti-apoptotic genes in MMterm·dom cells. In the unstimulated state, NF-κBterm·dom heterodimers are always kept in the IκBα|NF-κBterm·dom complex in the cytoplasm. Upon the stimulation of TNF, IKKterm·dom is transformed from its neutral state (denoted as IKKn) into its active state (denoted as IKKa). IKKa can phosphorylate and ubiquitinate IκBα and then release free NF-κBterm·dom heterodimers from the IκBα|NF-κBterm·dom complex. The free NF-κBterm·dom heterodimers enter the nucleus and regulate the transcription of IκBα, A20 and GADD45βterm·dom [ 24 ]. In addition to IκBα, A20 also serves as an inhibitor of the NF-κBterm·dom pathway, promoting the transformation of IKKa into its unactivated IKKi form [ 25 ]. As a downstream factor of NF-κBterm·dom signalling, the newly produced GADD45βterm·dom is expressed at a high level in the MMterm·dom cells and interacts with MKK7term·dom in the JNKterm·dom pathway, promoting the survival of MMterm·dom cells [ 4 , 26 , 27 ] (detailed information is included in part B). It should be noted that IκBβ, IκBγ and IκBε can also bind to NF-κBterm·dom heterodimers and inhibit their nuclear localization signals [ 28 ]. Only IκBα is considered in our model, as NF-κBterm·dom heterodimers are primarily sequestered by IκBα.

Part B describes how the overexpression of GADD45βterm·dom acts with MKK7term·dom and dysregulates the JNKterm·dom pathway to facilitate the survival of MMterm·dom cells. There are three cascades, JNKterm·dom (a member of the mitogen-activated protein kinase (MAPK) family), MKK7term·dom (a member of the MAPK kinase family) and MAPK kinase kinase (MAPKKK) in the JNKterm·dom pathway [ 8 , 29 , 30 ]. The activated kinase at the upper level phosphorylates the kinase at the next level down the cascade. TNF stimulation promotes the phosphorylation of MAPKKK mediated by reactive oxygen species (ROS) [ 26 , 31 , 32 ]. MAPKKK_P denotes phosphorylated MAPKKK. The full activation of both MKK7term·dom and JNKterm·dom requires the phosphorylation of two conserved sites [ 33 ]. MAPKKK_P promotes the single phosphorylation of MKK7term·dom (denoted as MKK7term·dom_P) at first and then doubly phosphorylated (denoted as MKK7term·dom_PP). Similarly, MKK7term·dom_PP then facilitates the single and double phosphorylation of JNKterm·dom (denoted as JNKterm·dom_P and JNKterm·dom_PP, respectively). JNKterm·dom_PP is able to promote the apoptosisterm·key of MMterm·dom cells by up-regulating pro-apoptotic genes [ 15 ]. The GADD45βterm·dom produced during NF-κBterm·dom signalling binds to MKK7term·dom and inhibits its enzymatic activity, which limits the phosphorylation of JNKterm·dom promoted by MKK7term·dom_PP. The decrease in JNKterm·dom_PP inhibits the expression of pro-apoptotic genes and MMterm·dom cell death.

Model equations To investigate the role of GADD45βterm·dom in MMterm·dom pathogenesis, the mathematical model proposed in [ 23 ] is extended by including the pathways of NF-κBterm·dom and JNKterm·dom as well as their crosstalk. The proposed model consists of 29 ordinary differential equationsterm·mth (ODEsterm·mth) with five state variables, namely osteoblast precursors (OBp), active osteoblasts and osteoclasts (OBa and OCa), MMterm·dom and the bone volumeterm·dom (BV). Equations (2.1)–(2.5) represent temporal variations in the concentrations of OBp, OBa, OCa, MMterm·dom and the BV, respectively. Equations (B1)–(B11) and (C1)–(C16) in the electronic supplementary material, describe variations in biochemical factors regarding the JNKterm·dom and NF-κBterm·dom pathways with time, respectively.

‘Hill functionsterm·mth’ are used to represent the cellular interaction via single ligand to receptor binding denoted by π functions [ 34 ], with equations (2.6) and (2.7) describing the stimulating and inhibiting functions of ligand-receptor binding. Here. L represents the concentration of the ligand, β represents maximal expression level of the promoter, n is the coefficient which regulates the steepness of the function π and k 1 and k 2 represent the dissociation constant, respectively. Both β and n are assumed to equal 1 in the model following the work of Pivonka et al . [ 34 ].

and dOB p /d t , dOB a /d t , dOC a /d t , (d/d t ) MMterm·dom and (d/d t ) BV denote the variations of OB p , OB a , OC a , MMterm·dom and BV, respectively. For example, dOB p /d t is the variation of OB p with time. D O B u and D O B p represent the differentiation rates of uncommitted osteoblast progenitors and osteoblast precursors. OB u and OB p are concentrations of uncommitted osteoblast progenitors and osteoblast precursors. π a c t , O B u T G F β represents the stimulation of uncommitted osteoblastic progenitors into osteoblastic precursors. π r e p , O B p T G F β represents the inhibition of the differentiation of osteoblastic precursors into active osteoblasts. π r e p , O B p V C A M 1 represents BMSC-MM cell adhesion that blocks the differentiation of mature osteoblasts from their progenitors. Definitions of other Hill functionsterm·mth and variables are not repeated here but they are included in electronic supplementary material, Appendix A.

The newly added Hill functionterm·mth π a c t , M M J N K _ P P in equation (2.4) represents the promotion of MMterm·dom cell apoptosisterm·key by JNKterm·dom_PP. The definition of π a c t , M M J N K _ P P is as follows: where K J N K _ P P represents the activation coefficient related to JNKterm·dom_PP promoting MMterm·dom cell apoptosisterm·key, and its value is included in table 1 . JNKterm·dom_PP represents the concentration of doubly phosphorylated JNKterm·dom. The calculation of JNKterm·dom_PP requires the mathematical modelling of the NF-κBterm·dom and JNKterm·dom pathways together with their interaction through GADD45βterm·dom. In accordance with the earlier work of [ 35 ] and [ 25 ], equations (B1)– (B11) and (C1)–(C16) in the electronic supplementary material, were constructed to simulate the JNKterm·dom and NF-κBterm·dom pathways. The distinct feature of these equations is that the interaction of the two pathways, which was ignored before, was included, with the addition of E MKK 7 GADD 45 β in electronic supplementary material, equations (B10) and (C16). E MKK 7 GADD 45 β represents the binding of GADD45βterm·dom to MKK7term·dom, which decreases the amount of MKK7term·dom available for phosphorylation into MKK7term·dom_P. This action further leads to the inhibition of the JNKterm·dom pathway and promotes the survival of MMterm·dom cells as discussed before. The definition of E MKK 7 GADD 45 β is as follows: where the first term represents the formation of the GADD45βterm·dom|MKK7term·dom complex, and the second term describes its dissociation. The definitions of the variables in the JNKterm·dom and NF-κBterm·dom pathways are included in electronic supplementary material, Appendices B and C. In summary, π function π a c t , M M J N K _ P P and 24 ODEsterm·mth are new additions to our model, allowing us to investigate the essential role of the interactions between NF-κBterm·dom and JNKterm·dom signalling in the development of MMterm·dom.

Descriptions and values of parameters used in the model.

Results The definitions and values of the model parameters are listed in table 1 . Several parameter values were reported in previous studies, while the remaining unknown parameters (i.e. those parameters where experimental data are unavailable or those that have no direct biological meaning) were fitted via the genetic algorithmterm·mth (GAterm·mth) in this paper. The initial values of the model variables are described in table 2 . The ode45 solver in the Matlab software package (R2014b, Mathworks, Natick, USA) is used to solve the model equations. Curves in figures 2 – 10 represent numerical solutions of model equations. The Matlab code is included in the electronic supplementary material. The calculation of unknown model parameters based on GAterm·mth are described as follows: and where X = [ K J N K _ P P , k G M , d G M , K D , V C A M 1 , M M , a c t , K D , S L R P s , M M , r e p , K D , I L 6 , M M , a c t , M M m a x ] is a row vector consisting of the seven unknown model parameters and represents one point in the parameter space. M ( X ) i and P i ( i = 1, 2, 3) represent model outputs corresponding to each point in the parameter space and the preferred model outputs, respectively. In this work, model outputs refer to the concentrations of OB p . Different combinations of parameter values have different outputs. GAterm·mth is an effective way to search for parameter values in parameter space to make model outputs approach to preferred outputs.

The initial values of cell concentrations used in the model.

Model simulations of the normalized variation in the concentrations of OBp, OBa, OCa and MMterm·dom cells with respect to their respective initial values (MMterm·dom cells are injected at 10th hour and removed at 150th hour).

Model simulation of the variation in NF-κBterm·dom due to the inhibition of binding of GADD45βterm·dom and MKK7term·dom, and the IKKterm·dom activation. ( a ) The binding of GADD45βterm·dom and MKK7term·dom inhibited by 40% at 150th hour. ( b ) The IKKterm·dom activation inhibited by 40% at 150th hour.

Figures 2 – 4 describe the variations in the concentrations of bone cells, MMterm·dom cells, the OBa : OCa ratioterm·dom, and the bone volumeterm·dom after the invasion and removal of the MMterm·dom cells. As demonstrated in figure 2 , the concentrations of OBp, OBa and OCa remain in a steady state under normal conditions. This steady state is disturbed due to the invasion of the MMterm·dom cells. OBp, OBa and OCa all increase to different degrees following the increase of MMterm·dom cells. These simulation results are consistent with experimental observations [ 44 – 47 ]. OCa grows in a larger degree than that of OBa, which leads the OBa : OCa ratioterm·dom to drop after the invasion of the MMterm·dom cells, as demonstrated in figure 3 . The OBa : OCa ratioterm·dom is a key factor in the variation of the bone volumeterm·dom, since OBa and OCa are for forming and resorbing bone, respectively, during the bone remodelling process. The decreased OBa : OCa ratioterm·dom causes the decline in the bone volumeterm·dom, as shown in figure 4 , which is consistent with the observed loss of the bone volumeterm·dom in MMterm·dom patients [ 44 ]. These simulation results match the work of [ 23 ]. Figures 5 and 6 reveal the change in the GADD45βterm·dom concentration after the activation of NF-κBterm·dom signalling and its influence on JNKterm·dom signalling. Figure 7 reveals the temporal variation in MMterm·dom cells due to the inhibition of GADD45βterm·dom and MKK7term·dom binding by 10% and 40%. Figures 7 and 8 suggest how GADD45β-targeted therapy inhibits the growth of MMterm·dom cells and JNKterm·dom_PP. Thus, the inhibition of GADD45βterm·dom binding to MKK7term·dom leads to an obvious drop of MMterm·dom cells, which then leads to a new steady state at a lower level. The larger degree of inhibition results in a larger drop in MMterm·dom cells. These simulation results can be confirmed by the experimental data [ 4 ]. Figure 9 shows the variation in MMterm·dom cells after IKKterm·dom activation is inhibited by 10% and 40%. As shown in figure 9 , IKKterm·dom inhibition leads to a decrease in MMterm·dom cells which is accentuated with a larger degree of IKKterm·dom is inhibition. Figure 10 indicates how the inhibition of GADD45βterm·dom and IKKterm·dom influence NF-κBterm·dom.

Model simulation of the variation in the normalized bone volumeterm·dom with respect to its initial value (MMterm·dom cells are injected at 10th hour and removed at 150th hour).

Model simulation of the variation in the normalized ratio of OBa to OCa with respect to the initial ratio (MMterm·dom cells are injected at 10th hour and removed at 150th hour).

Model simulation of the variation in the GADD45βterm·dom concentration after the activation of the NF-κBterm·dom pathway at 10th hour.

Model simulations of the variation in biochemical factors involved in JNKterm·dom signalling after the pathway is activated at 10th hour.

Model simulations of the variation in MMterm·dom cells due to the inhibition of GADD45βterm·dom and MKK7term·dom binding by 10% and 40% at 150th hour.

Model simulations of the variation in JNKterm·dom_PP due to the inhibition of GADD45βterm·dom and MKK7term·dom binding by 10% and 40% at 150th hour.

Model simulations of the variation in MMterm·dom cells due to the inhibition of IKKterm·dom activation by 10% and 40% at 150th hour.

Discussion The proposed model in this paper aims to investigate how the interaction between NF-κBterm·dom and JNKterm·dom signalling via the GADD45βterm·dom|MKK7term·dom complex promotes the development of MMterm·dom. The model can not only simulate the temporal variation of bone cells, the bone volumeterm·dom and biochemical factors involved in the NF-κBterm·dom and JNKterm·dom pathways, but it can also mimic the underlying mechanism in which the coupling between the two pathways inhibits apoptosisterm·key in MMterm·dom cells. Additionally, the efficacies of GADD45β- and NF-κB-targeted therapies were compared based on the model simulation. This work only demonstrates the temporal variation in the GADD45βterm·dom concentration after the activation of the NF-κBterm·dom pathway. The variations in other biochemical factors related to NF-κBterm·dom signalling were not repeated here as they have been analysed in [ 25 ] but are included in electronic supplementary material, Appendix C.

As mentioned above, JNKterm·dom_PP plays an essential role in the apoptosisterm·key of MMterm·dom cells, and JNKterm·dom_PP is able to inhibit the growth of MMterm·dom cells. According to figure 8 , the decrease in MMterm·dom cells can be explained by the rising amount of JNKterm·dom_PP after the inhibition of GADD45βterm·dom binding to MKK7term·dom. NF-κB-targeted therapies are most often performed by inhibiting IKKterm·dom, which is essential for NF-κBterm·dom activation [ 8 ]. Based on model simulations ( figure 10 ), the GADD45β-targeted therapy is suggested to be more sensitive than the NF-κB-targeted therapy as it can produce a greater suppression of MMterm·dom cells than IKKterm·dom inhibition when each are suppressed to the same degree. In addition to the higher sensitivity, GADD45β-targeted therapy has a high MMterm·dom cell specificity. This approach can not only kill MMterm·dom cells effectively, but it also has no effect on NF-κBterm·dom levels. Therefore, it can be predicted that the side effects seen with NF-κB-targeted treatment would be avoided. However, these in silico hypotheses generated from the mathematical model described in this paper will ultimately have to be validated by further experimental data to confirm the clinical potential of GADD45β-targeting therapies. It should also be noted that the proposed simplified model does not take into account the interactions with several other proteins and other cell structures, which represents a limitation of this study. Specifically, although the differentiation of progenitors into active osteoclasts and osteoblasts contains several intermediate stages, the model only considered four osteoblastic and three osteoclastic lineages, and contained three state variables: osteoblast precursors, active osteoblasts and active osteoclasts. Additionally, IκBs, including IκBα, IκBβ, IκBγ and IκBε, can also bind to NF-κBterm·dom and form a complex but we only considered IκBα.

Conclusion This paper describes a mathematical model of MMterm·dom that was used to investigate how NF-κBterm·dom signalling acts to include the interplay between NF-κBterm·dom and JNKterm·dom signalling that was not included in the previous models. The model not only reconstructs how the invasion of MMterm·dom cells disturbs the steady state of the bone microenvironment and triggers the variation in bone cells, but it also mimics the changes in biochemical factors involved in the NF-κBterm·dom and JNKterm·dom pathways. In addition, the model also investigates the efficacies of GADD45β- and NF-κB-targeted treatments, suggesting that GADD45β-targeted therapy is more effective but has no apparent side effects. The simulation results match the experimental observations. This model helps to illuminate the essential function of the crosstalk between the NF-κBterm·dom and JNKterm·dom pathways during MMterm·dom development. It is anticipated that this model could be employed as a useful tool to initially investigate and even explore potential therapies involving the NF-κBterm·dom and JNKterm·dom pathways. More work is required to improve the built model to take into account many others parameters such as the interactions with other protein partners, with membranes and other cell structures.

Supplementary Material Supplementary Material Supplementary Material Supplementary Material Supplementary Material Supplementary Material Supplementary Material Supplementary Material Supplementary Material Data accessibility Matlab code is available in the electronic supplementary material.

Authors' contributions Y.Z. carried out program coding and model simulations, participated in the design of the study and drafted the manuscript; C.Z. carried out the analysis of experimental data and critically revised the manuscript; Q.Y. collected the experimental data and critically revised the manuscript; B.J. conceived of the study, designed the study, coordinated the study and helped draft the manuscript. All authors gave final approval for publication and agree to be held accountable for the work performed therein.

Competing interests We declare that we have no competing interests.

Funding This work was supported by National Key R&D Program of China (grant no. 2018YFB1305400); the National Natural Science Foundation of China (grant nos. 61673246 and 81301294); the key Research and Development Program of Shandong province (grant no. 2016GSF201168) and the Research and Development Program of Jinan (grant no. 201907064).

References Associated Data Supplementary Materials Data Availability Statement Matlab code is available in the electronic supplementary material.