CC BY
16
1 Institute of Computational Technologies of SB RAS 630090, Novosibirsk, Lavrent'ev ave., 6
2 Institute of Molecular Pathology and Pathomorphology 630117, Novosibirsk, Timakov str., 2
This paper is devoted to a numerical analysis of the solutions of the equations system describing the dynamics of the concentration of p53 and Mdm2 proteins in their interaction. The detailed study of the solutions was made when the mathematical model parameters deviated from the basal values. The system states in which the threat of the cell uncontrolled apoptotic death accelerated the organisms aging processes or the excessive suppression of p53-induced apoptosis increased the tumors risk have been investigated within the numerical experiments. The mechanism of the system managing under stress conditions has been studied.
Keywords: tumor markers; p53-Mdm2-network; time delay equations; numerical modelling.
удк 519.6; 519.8; 51-7; 616
MATHEMATICAL MODELLING OF THE TUMOR MARKERS NETWORK Olga Falaleevna VOROPAEVA1, Yurij Ivanovich SHOKIN1,
Lev Moiseevich NEPOMNYASHCHIKH2 , Svetlana Robertovna SENCHUKOVA2
The p53 protein (tumor necrosis factor), involved in many life and death processes, including the formation of tumors and aging, is expressed in all the cells of the organism. Mdm2 protein is considered to be the key negative p53 regulator. The mechanism for the functioning of the p53-Mdm2 system is very complex, so the answers to questions that are important in clinical practice are often found to be mutually exclusive. In particular, it has been found that the imbalance in the p53 and Mdm2 interaction can be the cause for serious pathological changes in organs and tissues. Thus, the consequences of the faults, caused by the excessive production of p53 and the superactivation of the p53-dependent apoptosis, are often serious diseases such as neurodegeneration (Alzheimer's, Parkinson's, multiple sclerosis, epilepsy), osteoporosis, growth arrest, and premature aging of the internal organs. At the same time, the loss of the p53 protein's function is found in approximately 50 % of human malignant tumors. Therefore, the inactivation of p53 is also conventionally seen as an undesired and dangerous event, while in antitumor therapy the restoration of an active p53, including as a result of the artificial break of the interaction between p53 and Mdm2, is thought of as one of the key elements ensuring the
normal passage of signals stopping growth stop and/or apoptosis in cancer cells. It is necessary also to take into account that p53-dependent apoptosis induced in normal tissues can become the cause for serious side effects, thus limiting the efficiency of the antitumor therapy. In this connection, the artificial inhibition of p53 in certain tumors, especially in those without functional p53, can have a protective effect, supporting the normal functioning of the tissue in the process of regeneration.
Thus, the investigation of the function of the p53 and Mdm2 proteins and of the mechanism of their interaction is paramount both for developing new approaches to cancer treatment and determining the prevention strategy for many diseases, including measures to slow the aging processes. Many aspects of this problem have so far remained little known and thus are now at the center of attention among researchers [1-3, 5] (these works contain detailed bibliographies). Under these conditions, the mathematical simulation of the interaction of the p53 and Mdm2 proteins is one of the effective and accessible methods of analysis of the development of the p53-Mdm2 system and of exits from peak situations.
Voropaeva O.F. - doctor of physical mathematical sciences, senior researcher, e-mail: vorop@ict.nsc.ru Shokin Yu.I. - academician of RAS, director
| Nepomnyaschchikh L.M. \ - corresponding member of RAS, director Senchukova S.R. - doctor of medical sciences, leading researcher
MATHEMATICAL MODELS AND NUMERICAL
TECHNOLOGIES
We consider in this work two interrelated mathematical models of the p53-Mdm2 network. The first (basic) model of the proteins concentrations dynamics includes the system of two nonlinear equations with the retarded argument [4] (see also [5, 6]):
dyi/dt= s - a f oi(t), y2(t), kf) - byl(t), (1) dy2/dt = cig(y}(t - t), y2(t - t), f, kg) - c^co, (2)
where the interaction of the proteins is determined by nonlinear functions f and g. Here y1 and y2 are concentrations of the p53 and Mdm2, respectively; s is the rate of p53 production; a is the rate of p53 degradation by means of ubiquitination; b is the rate of p53 spontaneous degradation; c1 is the rate of Mdm2 production, including that by its interaction with p53; c2 is the rate of Mdm2 protein degradation; kf is the constant of the p53-Mdm2 complex dissociation constant; and kg is the p53 protein and Mdm2 gene dissociation constant (a more detailed biological description of the parameters is found in [4]). The initial data for system (i)-(2) are given as the functions of «history» ^k(9): yk(9) = ^k(9), 0 e [-t, 0], k = 1,2. The basal values of the parameters [4] are coordinated with the data from different laboratory studies. The time delay parameter t determines in the Mdm2 reaction to a change in the state of the p53 protein (more detail see in [1-3]).
The second model describes hypothetical stages of process and uses the simplest ODE system of sufficiently higher dimension. We show numerically that in the passage to the limit in which the second model has infinitely many stages we obtain model based equation with retarded argument.
The characteristic features of the present model are nonlinearity, time delay, and the possible rigidity of the system and bifurcation of solution under some combinations of parameters. A sufficiently detailed description of the algorithm's testing procedure is present in [6]. The solution of problem (1)-(2) is based on the widely known method of steps. From considerations of the simplicity of the numerical implementation, in all the calculations the delay value t was taken as a multiple of the step of the computational grid (on condition that the grid step is a rather small value, this limitation avoids having a substantial effect on the character of the solution). At each time step, the problem's solution is found by numerical methods with the attraction of iterations by nonlinearity. For the sake of convenience of the numerical implementation, the system of equations (1)-(2) can be nondimensionalized with the use of the delay value t as the time scale. This procedure is equivalent to compressing the time
axis t times. In doing this, one should take into account that in transition to a fairly high t the accumulation of errors caused apparently by rigidity is typical for the system in the transformed form. To exclude undesirable computational effects, the choice of the numerical method and controlling the influence of the grid step on the solution have been given the closest attention. The methodical calculations of system (1)-(2) were performed on a sequence of grids with a resolution ranging from 1 to 10000 points per second. Based on these tests, in the numerical experiments, the grid with a resolution of not more than 100 points per second was acknowledged as optimal.
Based on the results of the test calculations and comparison of a sufficiently high number of known numerical methods for the solution of problem, the explicit Adams method with iterations by nonlinear-ity was used. In doing this, the defect of the multistep methods related to the need for a special approach to the simulation of the solution at the first points, in case of equations with delay, is compensated by calculations on sufficiently high time intervals.
In [4-6] it is shown that under basal conditions system (1)-(2) has stationary solution only as im-
t t t • •, • , / basal basaU , 1 n
movable limit point (y1 , y2 ) at any values of delay, including t = 0, and in the disturbance of the basal state points the Andronov-Hopf bifurcations can be observed. All the approximated coordinates of the immovable points are with sufficiently high precision in keeping with the analytical stationary solution [4].
Attention was given to the question of the influence of the choice of the initial data on the solution of the problem (1)-(2). A series of numerical experiments were made, in which either zero values or constant stationary basal values of concentrations y1 and y2 with the introduction of stress perturbations in them ranging from 10 % to 4 orders (both decreasing and increasing, alternately for each component of solution y1 and y2 or simultaneously for both) were used as the «history» functions ^k(9). In doing this, the time delay parameter was also varied, and the values of the parameters of the system s, a, b, c1, c2, kf and kg were taken as basal. In addition, the numerical experiments, in which as the initial data on interval [-t, 0] the periodic solutions of this problem were used, and hereinafter at all the calculations were made under basal conditions. The calculations show that in all t > 0 the considered stress situations caused by the change of the initial level of protein concentration, the character of the solution of system (1)-(2) does not practically change independent of the value of the delay and the initial concentration levels: y^t) and y2(t) returns to the earlier found limit yj3sa, y2asa; i.e., the solution of the problem has a stable focus. In using as the initial data the
, , • t n , basal t basal
stationary values of concentrations y and y2 with the introduction in them of stress perturbations ranging from 20 % to 2 orders (towards a decrease and increase, alternately for each component and simultaneously for both, the delay parameters also varied) similar results were obtained.
In addition, numerical experiments were performed, in which different numerical methods were used for the reproduction of a periodic solution (delay t reached the bifurcation value tbif). The stationary solution is the focus corresponding to a rather small t < tbif was used as the initial data. The comparison of the solutions obtained with the use of the predictor-corrector method with iterations by nonlinearity and the Runge-Kutta method of order 4 has shown that both numerical methods yield rather close solutions, and the introduction of the iterations allows a considerable increase in the accuracy of the calculation by using of the predictor-corrector method. The solutions obtained by use of the Adams method of order 4 and those by the Runge-Kutta method of the same order were found to be the closest. In the calculations, where the obtained periodic solutions were used as the initial data, a comparison of the methods yields similar results. Then the real order of accuracy of the method on such solutions ranges from 1 to 2.83.
DEPENDENCE OF THE SOLUTIONS ON THE
CHANGE OF THE MODEL PARAMETERS
We have carried out a numerical study of the reaction of system to stresses in the form of a deviation of the parameters from the basal values. The numerical experiments were made in a sufficiently large interval of the delay values ranging from a few seconds to many hours, which is characteristic of the given biological system (see, for example, [1-3]).
Fig. 1 provides some data on the obtained solutions that correspond to the alternatively deviation of the parameters from the basal values. It has been found that the basic part of the solutions are immovable points (focuses). Moreover, the calculated coordinates of the immovable points in the phase plane agree, with a rather high degree of accuracy, with the analytical stationary solution of problem [4], i.e., with the solution of system (1)-(2) at the corresponding values 5, a, b, c1, c2, kf and kg. In the course of these numerical experiments, it has been found that the dissociation constants kf and kg become bifurcation parameters, providing the start of periodic oscillations in the p53-Mdm2 network, at the time of delay typical of real conditions.
The analysis of the numerical data shows that within the accepted model, in the deviation of one of the parameters of the system from the basal value, there may appear the following typical situa-
tions: a) a considerable accumulation of p53 at a rather low level of its inhibitor Mdm2; b) a considerable accumulation of Mdm2 at a lower level of p53; c) rapid attenuation of the oscillations of protein concentrations with the achievement of constant stationary values close to the basal ones; and d) periodic oscillations. The calculations show that time delay is one of the key parameters affecting the state of the p53-Mdm2 system. The attenuating solutions are typical for most of the arising situations, which testifies to the stable reaction of the system to the deviation of the parameters from the normal values. We note that although the periodic oscillations of the concentrations in the p53-Mdm2 system, with regard to some characteristic constant values, which are often rather close to the basal values; however, the maximum values of the concentrations achieved in this may become critical for the normal functioning of the system and the organism as a whole.
Note that the obtained data do not conflict with the known ideas on the phenomenon being simulated and allow the conclusion on the adequacy of the considered mathematical model of the p53-Mdm2 grid. In particular, the presence of oscillations of the p53 and Mdm2 concentrations taking place almost in the counterphase are in correspondence with the known data of the laboratory experiments (see, for example, [1-3] and references in [5]).
RESTORATION OF THE P53-MDM2 FEEDBACK
The analysis of the data presented in the previous section shows that a considerable part of the
Fig. 1. The stationary solutions of the problem with the alternating variation of parameters: solid lines -immovable points (1, 1', 2) and the limit cycles (3, 4) of the numerical solutions; dashed lines -analytical stationary solutions [2]; 5 - the stationary solution under basal conditions
solutions of system (1)-(2) (including solutions with an immovable point) correspond to situations in which the disturbance of the feedback of proteins p53 and Mdm2 is characteristic. The negative and positive consequences of these disturbances in real medico-biological processes have been discussed above. In this section we will consider in greater detail the most characteristic of the noted situations and present the results of the numerical experiments aimed at finding solutions that are more favorable for the live organism-understood as the achievement of the basal values of the levels of concentration of proteins.
The methodology of the numerical experiments in this section is as follows: we will try to compensated disturbances of the feedback in the p53-Mdm2 network that are the result of a stress for some parameter of system (1)-(2) by the repeated stress for one of the remaining parameters. This methodological technique is fairly well known in medico-biological studies as suppression, i.e., total or partial suppression of the first (direct) mutation by the second (suppressor) mutation, which may lead to the partial or total restoration of the normal phenotype. Our purpose is to find that parameter of the mathematical model, whose change will lead to the restoration of the basal state of the p53-Mdm2 system. As the considered states of the system are characterized by solutions with an immovable limit point, for which earlier stability relative to the initial
Fig. 2. Problem solutions at a = 0.5abasal: 1 - basal state, 2 - s = 0.5sbasal, 3 - b = 35bbasal, 4 - kf = = 0.13kfbasal (t = 120 s); 5 - kf = 0.15kbasal (fragment of the curve in the vicinity of immovable point), 6 - limit cycle at kf = 0.13kbasal and t = = 25000 s; the dashed lines show the analytical stationary solution [5] at basal values of s, b, kf, kg, c1, c2
data was found, for simplicity of the calculation procedure, all the numerical experiments are made at zero initial conditions. All the calculations have been continued down to the delay values of 30 000 s and the duration of the observed processes were not less than 120 days.
Excessive p53 level
Within the considered mathematical model, the excessive p53 level threatening excessive p53 dependent apoptosis can be caused by the lack of balance in the rate of p53 generation and Mdm2 degradation, as well as by disturbances in the mechanism of the interaction of proteins, which is regulated through the introduction of the dissociation constants.
For example, a considerable excess in the value of the rate of p53 generation (parameter s) leading to threatening increase of the p53 level can be, according to the calculations, compensated by increasing the rate of the p53 degradation through ubi-quination a or of the rate of p53 spontaneous decomposition b, as well as by a lower dissociation constant of the p53-Mdm2 complex kf. When a decreases more than by one order, the situation is complicated with a further decrease of p53, and the upper boundary of the Mdm2 concentration is very close to the basal value. The increase of kf also leads to the further growth of p53. The decrease of b will not change the situation. The decline of the rate of the p53 degradation through ubiquitination (parameter a), which is able to cause the disturbance of the feedback of p53 and Mdm2, can be compensated within the considered model mostly due to the variation of the parameters s, b, kf (see fig. 2). Variation of the parameters ch c2, and kg can compensate the negative situation only partially, since in this case the levels of protein concentrations exceeding the basal values by a factor of at least 1.5 to 2 are installed.
Thus, in order to remove the threat of the uncontrolled apoptotic death of cells, first of all it is possible to apply the connection of two pairs of parameters - the rate of p53 generation and the rate of this protein degradation through ubiquitination or the Mdm2 degradation rate and the rate of its generation: an undesirable increase or decrease of one of the parameters in each of these pairs is compensated by a similar (by the same factor) change of the second parameter of this pair. In this situation, the system is considerably less sensitive to both the regulation through the rate of the spontaneous p53 decomposition and also to the strengthening of the p53 and Mdm2 bond through dissociation constants. At sufficiently large values of the t delay (close to the observed ones), it is not easy to remove the
negative consequences caused by the increased rate of p53 generation by the control of the dissociation constants.
Excessive Mdm2 level
In the course of the numerical experiments, there was a situation in which there was an increased risk of a decrease of the p53 function, when the interaction of p53 and its inhibitor Mdm2 is attenuated, and as a result the p53 level is found to be very low and the level of Mdm2 itself is too high. Within the considered mathematical model, this state arises when the rate of Mdm2 degradation is low or the rate of its generation is too high. We are considered the ways to compensate these disturbances.
The damping of the interaction in the p53-Mdm2 network, accompanied by the fast growth of the Mdm2 concentration, can be compensated within the considered model only through controlling the rates of Mdm2 generation or degradation, and at very moderate delay values it can be also done by a significant increase of the dissociation constant of the p53 protein and the Mdm2 gene. The studies described in this section allow us to answer important questions for the clinical practice on the functioning and regulation of the p53-Mdm2 grid, including those associated with the artificial break of the p53-Mdm2-interaction for therapeutic purposes and with the actions for the restoration of the basal state of the system.
CONCLUSIONS
Balancing the rates of cellular proliferation and apoptosis is an important tool for maintaining the homeostasis of the normal tissue. Disturbance of this balance can become a cause of serious pathological changes in organs and tissues. In this work, a numerical investigation of the disturbances in the system of p53 and Mdm2 proteins, which participate in the correction of DNA defects and are involved in many vital processes, including onco-genesis, neurodegenerative diseases, and aging, was carried out. We are considered two interrelated mathematical models of the p53-Mdm2 network. The basic model of the proteins concentrations dynamics includes the system of two nonlinear equations with the retarded argument. The second model describes
hypothetical stages of process and uses the simplest ODE of higher dimension. We show numerically that in the passage to the limit in which the second model has infinitely many stages we obtain model based equation with retarded argument.
All obtained numerical solutions have a sufficiently clear biomedical meaning. The adopted mathematical models describe the functioning of the p53-Mdm2 network in normal conditions and predict possible dangerous situations for the organism. Stress situations associated with the emergence of an imbalance in the rates of p53 and Mdm2 generation and degradation and also with disturbances in the mechanism of interaction of proteins regulated within the considered models through the constants of dissociation have been studied. The variants of the compensation of the considered disturbances by repeated stress for the system parameters have been numerically investigated.
ACKNOWLEDGEMENTS
The work was performed with the financial support of the Council on Grants (grant NSh 5006.2014.9) of the President of the Russian Federation.
REFERENCES
1. Batchelor E., Loewer A., Mock C., Lahav G. Stimulus-dependent dynamics of p53 in single cells // Mol. Syst. Biol. 2011. 7. ID 488.
2. Lane D., Levine A. p53 Research: the past
thirty years and the next thirty years // Cold Spring Harb. Perspect. Biol. 2010. 2. (12). ID a000893.
3. Loewer A., Batchelor E., Gaglia G., Lahav G. Basal dynamics of p53 reveal transcriptionally attenuated pulses in cycling cells // Cell. 2010. 142. (1). 89100.
4. Tiana G., Jensen M., Sneppen K. Time delay as a key to apoptosis induction in the p53 network // Eur. Phys. J. B. 2002. 29. (1). 135-140.
5. Voropaeva O.F., Shokin Yu.I., Nepomnyash-chikh L.M., Senchukova S.R. Mathematical model of function and regulation of p53-Mdm2 biological system. Moscow: Publishing House RAMS, 2014 [In Russian].
6. Voropaeva O.F., SenchukovaS.R., BrodtK.V. et al. Numerical simulation of ultradian oscillations in p53-Mdm2 network under stress conditions // Math. Models Comput. Simul. 2015. 7. (3). 281-293.
удк 519.6; 519.8; 51-7; 616
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ СЕТИ ОНКОМАРКЕРОВ
Ольга Фалалеевна ВОРОПАЕВА1, Юрий Иванович ШОКИН1,
Лев Моисеевич НЕПОМНЯЩИХ2, Светлана Робертовна СЕНЧУКОВА2
1 ФГБУНИнститут вычислительных технологий СО РАН 630090, г. Новосибирск, пр. Академика Лаврентьева, 6
2 ФГБНУ Институт молекулярной патологии и патоморфологии 630117, г. Новосибирск, ул. Тимакова, 2
В работе выполнен численный анализ решений системы уравнений, описывающих динамику концентраций белков р53 и при их взаимодействии. Проведено детальное численное исследование решений при
отклонении параметров математической модели от базальных значений. В рамках численных экспериментов получены состояния сети p53-Mdm2, соответствующие как угрозе неуправляемой апоптотической гибели клеток, ускоряющей процессы старения организма, так и ситуации чрезмерного подавления апоптоза, когда увеличивается риск онкологических заболеваний. Исследованы механизмы управления системой p53-Mdm2 в условиях стрессов.
Ключевые слова: онкомаркеры, сеть p53-Mdm2, уравнения с запаздывающим аргументом, численное моделирование.
Воропаева О.Ф. - д. ф.-м.н., старший научный сотрудник, e-mail: vorop@ict.nsc.ru Шокин Ю.И. - академик РАН, директор | Непомнящих Л.М. | - член-кор. РАН, директор Сенчукова С.Р. - д.м.н., ведущий научный сотрудник
