A new treatment for breast cancer using a combination of two drugs: AZD9496 and palbociclib

In this study, we examined a mathematical model of breast cancer (BC) treatment that combines an oral oestrogen receptor inhibitor, AZD9496 with Palbociclib, a selective inhibitor of cyclin- dependent kinases CDK4 and CDK6. Treatment is described by analytical functions that enable us to control the dosage and time interval of the treatment, thus personalising the treatment for each patient. Initially, we investigated the effect of each treatment separately, and finally, we investigated the combination of both treatments. By applying numerical simulations, we confirmed that the combination of AZD9496 with palbociclib was the optimal treatment for BC. The dosage of AZD9496 increased and decreased throughout the treatment period, while the intervals were constant between treatments. Palbociclib changed almost cyclically, whereas the time intervals remained constant. To investigate the mathematical model, we applied the singularly perturbed homotopy analysis method, which is a numerical algorithm. The significant advantage of this method is that the mathematical model does not have to contain a small parameter (as is standard in perturbation theory). However, it is possible to artificially introduce a small parameter into the system of equations, making it possible to study the model using asymptotic methods.

The growing awareness of the disease, responsiveness to early detection, and improvements in planning and carrying out treatments according to the characteristics of the tumour have led to a significant increase in the cure rate.
To investigate the nonlinear ODE system, we combined two semi-analytical methods: the homotopy analysis method (HAM) and the method of integral invariant manifold (MIM).HAM employs the concept of homotopy from topology to generate a convergent series solution for nonlinear systems.This is a mathematical technique for solving a system of differential equations without a small parameter appearing in the system 22 .
MIM is a technique used to simplify the study of reaction mechanisms using dynamic systems, as proposed in 23,24 .

Mathematical models of BC
In recent years, global scientists from various fields, such as physics and mathematics, have integrated their research into applications in medicine and cancer research in particular [25][26][27][28] .The main aim is to devise a mathematical model (usually a system of nonlinear ODEs or PDEs) that considers various parameters related to the human body, such as different cells from the immune system and the interaction of these parameters with the tumour itself and with drug treatment.By solving a mathematical model numerically or analytically, it is possible to obtain the size of the tumour as a function of time, which is one of the most significant pieces of data in cancer research.As can be seen from many diverse studies, mathematical models can successfully and accurately predict the size of tumours [29][30][31][32] .
There are several significant advantages of using mathematical models in cancer research.The main advantage is that if one wants to adjust the treatment personally (such as the dosage of the drug and the time of administration of the treatment), there is no need to perform experiments on animals or humans, but simply change the parameters of the mathematical model.Another advantage is that an increasing number of equations related to the aforementioned interactions can be added, making the mathematical model more accurate.For example, drug treatment can be described by an explicit function that controls the doses of the drugs and changes the treatment administration times.Another significant advantage is the combination of machine learning (ML) and mathematical models; therefore, a very large dataset can be considered.
Here, we present important examples of BC research using mathematical models.BC is the most common malignancy in western countries.When the disease is diagnosed in its early stages, the chance of a cure increases to approximately 90% .It appears that the number of women who recover from BC is constantly increasing owing to early detection, improvements in treatment methods, and the widespread increase in global awareness.An example of this success is a 2018 study that included 699 women aged 18-85 with AH; the researchers 33 created a three-variable model aimed at predicting the risk of developing BC among women with atypical hyperplasia (AH) at biopsy.The three variables were age at the time of biopsy, age at the time of biopsy squared, and the number of AH foci.Women with AH identified on breast biopsy have an increased cumulative risk of developing BC, but the risk prediction models that exist today do not provide an accurate personalised assessment of the risk in this group.In this study, we analysed breast disease (benign growth) to develop and validate a BC risk prediction model that is specifically tailored to women with AH, called AH-BC.Researchers have created a new model that could more accurately predict the risk of BC in women with AH.After 10 years, the AH-BC model demonstrated good discrimination and calibration and was also valid when tested on an external model.
The next example is a study presented by 10 which describes the combination of a mathematical model describing chemotherapy treatment for BC with a ML algorithm to increase the performance in predicting tumour size using a five-step procedure.The focus of this research involves extracting the numerical solutions obtained from the mathematical model and integrating them into ML as additional data.The authors of this study established this procedure in four cohorts of women of different ages with BC who received chemotherapy.A performance comparison among the algorithms showed that the root mean square error of the linear regression decreased with the addition of the mathematical results, and the accuracy of prediction and the F1-scores increased with the addition of the mathematical model to the neural network.
Another study on mathematical models for BC research was presented by 34 in 2021.The authors developed a mathematical model to investigate the interaction between immune cells and BC tumour development.The main conclusion of this study, based on the application of numerical simulations to a mathematical model, demonstrated the direct effects of macrophages and adipocytes on cancer cell growth.In addition, interactions among immune cells can affect the immunological responses in various cancer types.

Formulation of the system of equations
In this section, we present a mathematical model 35 for ER-positive BC treatment using two different treatments: AZD9496 and palbociclib.We propose a new personalised treatment based on analytical functions which depend on two parameters: the dosage of medicine and the time interval between treatments.These two parameters enable us to control the treatments such that we can change the dosage and time interval depending on the tumour size.For this purpose, we developed an ODE equation for the function of treatments as a function of tumour size.The solution profiles of these equations show the dosage and time interval as a function of the tumour size at each given time.
The dynamical variables of the model are as follows: C C [cell] is the MCF-7 tumour cell population, N K [cellL −1 ] , is the natural killer (NK) cell population, W BC [cellL −1 ] is the white blood cell (WBC) population, C TL [cellL −1 ] is the cytotoxic T cell (CTL) population, A nc ZD [mg] is the AZD9496 not in circulation, A c ZD [mg] is the AZD9496 in circulation, P nc a [mg] is the palbociclib not in circulation, P c a [mg] is the palbociclib in circulation, and F [mg] and H [mg] are functions of the treatment of AZD9496 and palbociclib, respectively ( q A ZD is the amount of AZD9496 and q P a is the amount of palbociclib).The mathematical model is a system of ordinary dif- ferential equations, nonlinear first order in the form of The initial conditions of the model at t = 0 are In the next section, we apply the singularly perturbed homotopy analysis method (SPHAM) to investigate the mathematical model above.

Combination of HAM and MIM
In this section, we present HAM, MIM, and a combination of these two semi-analytical methods.

HAM
In this section we introduce the HAM procedure as present in 36 .
Given a nonlinear system of differential equations in the form of where N is a nonlinear operator, r is a vector of spatial variables, t is time, and u is an unknown function.By generalising the traditional concept of homotopy 37 , constructed the so-called zero-order deformation equation where ℏ is a nonzero auxiliary parameter, H is an auxiliary function, ℓ is an auxiliary linear operator, u 0 (•) is the initial guess for u(•) , and is an unknown function.The degrees of freedom are used to choose the initial guess, auxiliary linear operator, auxiliary parameter, and auxiliary function, H.To obtain the corresponding mth-order deformation of the HAM, we define the vector as follows: (1)  14) m times with respect to the embedding parameter p, setting p = 0 , and finally dividing the terms by m!, we obtain the mth-order deformation equation in the form where, Applying the inverse operator ℓ −1 (•) to both sides of Eq. ( 16), we obtain Thus, it is easy to obtain u m for m ≥ 1 at mth-order and finally obtain the solution of ( 13) as follows: In our model, we choose the initial estimates to be θ g (0) = 0 and r m (0) = 1 which satisfy the initial conditions.The linear operator is as follows: with the property ℓ(c 1 τ + c 2 ) = 0 , where c 1 and c 2 are the constants of integration.Following the above procedure, we define the series' for the dynamic variables C C , N K , W BC , C TL , A nc ZD , A c ZD , P nc a , and P c a as Based on the above notation, we let The set of equations ( 21) defines the analytical solutions as a power series of the dynamic variables of the system (1)-( 8) and the set of equations ( 21) defines the linear operator.

MIM
In this section, we introduce the primary theory underlying MIM.To apply the MIM, the system of differential equations must be in the form of a singular perturbed system: where 0 < ǫ << 1 is the small parameter of the system, � x = (x 1 , . . ., x k fast ) ∈ R k fast , and � y = (y 1 , . . ., y l slow ) ∈ R l slow are the fast and slow dynamic variables, respectively, of the system ( k fast + l slow = n ).The functions T and Q are assumed to be C ∞ of x , y and ǫ in (V × U) × I where V and U are open subsets of R k fast and R l slow , respectively.
I is an open interval containing the small parameter ǫ.
In general, when dealing with models describing natural phenomena, the small parameter does not appear explicitly in the system of equations, but in a hidden form; that is, the hierarchy of the mathematical model is hidden.
To apply MIM, the equations for m deformations of the linear operator can be rewritten as follows: where F 0 , F 1 , and F 2 are the corresponding terms in Eq. (23).
We obtain the second order by the small parameter γ perturbation of the equation According to the standard regular perturbation theory, the previous equation is a zero-order approximation of the perturbed equation 37 (i.e., terms with γ and γ 2 vanish).The first-order approximation is as follows: This procedure enabled us to expand the fast variable to a power series of small system parameters.

Combination of HAM with MIM
In this section, we present a combination of a HAM and MIM called SPHAM.A combination of these two methods is necessary because the mathematical model presented in the system of Eqs.(1)- (11) does not have an explicit hierarchy; that is, the hierarchy of the system is hidden; therefore, we cannot apply MIM separately.However, HAM does not require a small parameter of the system; the small parameter artificially enters the system of differential equations, which is problematic because not all variables of the mathematical model are fast.Therefore, we developed an algorithm that combined these methods.
The algorithm for the combination of HAM with MIM: 1: Substitute γ = 0 = γ 2 into Eq. ( 24) 2: Substitute m = 1 in Eq. ( 17) (the expression of R m u , where u is the dynamic variable of the system).From the expression R m u 1 (where 1 is the first variable of the system), we obtain the expression u 0 2 (where 2 is the second variable of the system) 3: Substitute the initial condition u 0 2 = u 2 (0) in step 2 and obtain an equation for u 0 1 : 4: Analytically solve the equation for u 0 1 .5: Substitute u 0 1 and the initial condition u 0 2 = u 2 (0) into the expressions R 1 u 1 and R 1 u 2 , 6: Repeat steps 1 − 4 for the remaining variables and initial conditions of system (1)-( 11), 7: Substitute the expressions from Step 6 into the following 1-order deformation: 8: Repeat steps 1-7 m times to obtain the m-order deformation of each variable of the mathematical model.9: Obtain the following explicit expressions: where u represents the variables in the system.10: Sum the expressions from step 9 and obtain the following semi-analytical approximation for the dynamic variables of the mathematical model: .

Results and discussion
In this section, we investigate the mathematical model (1-11) with the initial conditions (12) which presents BC treatment combining an oral oestrogen receptor inhibitor called AZD9496 with palbociclib, which is a selective inhibitor of the cycling-dependent kinases CDK4 and CDK6 for non-standard protocols, using the functions F and H for different doses and time intervals.palbociclib was administered orally and was absorbed slowly from the intestine in 6-12 h.We applied a combination of HAM and MIM and obtained the semi-analytical solution profiles of the dynamic variables of the mathematical model ( 1)-( 11) using the parameters presented in Tables 1, 2.
Figure 1 shows the solution profiles of the cancer cells AZD9496 and palbociclib.AZD9496 is administered at a constant dose of 300 mg.One tablet daily for 14 days, followed by 14 days off, for four cycles of administration.Palbociclib is administered arbitrarily, i.e., the dose of the drug changes arbitrarily, and the time interval of the drug is constant.The maximum and minimum doses were 125 mg and 20 mg, respectively.In this case, as can be seen from the solution profile of the cancer cells, the cancer decreased from the initial condition to zero after 72 days.
Figure 2 presents the solution profiles of the cancer cells and drug protocols.The AZD9496 administered did not follow the standard protocol; that is, the dose of the drug was changed, whereas the time intervals were constant.The administration of AZD9496 is as follows: 300 mg every day for 14 days, break for 14 days, 135 mg every day for 14 days, break for 14 days, 250 mg every day for 14 days, break for 14 days, and 50 mg every day for 14 days.
Palbociclib is not a standard protocol.In this case, the drug dosage changed almost cyclically.The time intervals were kept constant.The dose of palbociclib starts at 125 mg, is reduced to 85 mg and then 75 mg, after which the dose again is increased to 125 mg, reduced to 100 mg and then 60 mg, again increased to 125 mg, and finally reduced to 110 mg.At the end of the treatment, the dose of palbociclib was reduced to 40 mg.
In this combination of treatments, we achieved optimal treatment for the cancer cells which were reduced to zero after 32 d of treatment, as can be seen from the solution profile of the cancer cells presented in Fig. 2.
Figure 3 shows the solution profiles of the cancer cells and drug protocols.The doses of AZD9496 and palbociclib were constant in terms of dosage and time intervals: 250 mg of AZD9496, a cycle of one tablet per day for 14 days and a break of 14 days.The dose of palbociclib was maintained at 125 mg at constant intervals.
In this case, the number of cancer cells decreased from the initial condition to zero after 72 days.The significant difference between Figs. 1 and 3 is that, although both decrease to zero after 72 days, in the last graph, the cancer cells decreases to zero more quickly than those in the previous graph.However, the dose of AZD9496 was altered.The dose did not change regularly; that is, it increased and decreased.At the beginning of the treatment, the doses of AZD9496 were 300, 225, 250, and 50 mg.The drug administration times were different from those in the previous cases.The drug was administered for 21 days, with a break of 7 days and then administered again for 21 days.
In this case, the solution profile of cancer cells decreased rapidly to zero after 72 days.The graph decreases approximately linearly to zero.
Figure 5 shows the profiles of untreated cancer cells.It is clear that in this case, cancer cells grow rapidly.Figure 6 presents the solution profiles of the cancer cells and drug protocols.The dose of palbociclib increased from 20 to 125 mg.However, the dose of AZD9496 constantly decreased, starting at 300 to 225, 125, and finally, 50 mg.In this case, we see that the cancer cells, at the beginning of the treatment, started to increase for 23 days from the initial conditions up to 1.2 • 10 8 and then started to decrease for 45 days until they reached zero.The biological explanation for this phenomenon is that at the beginning of treatment, the dose of palbociclib was very low (20 mg) and then gradually increased.
Figure 7 presents the solution profiles of the cancer cells and drug protocols.In contrast to Fig. 6, in this case, the dose of palbociclib constantly decreased, starting from 125 mg, whereas the dose of AZD9496 constantly increased from 50 to 125, 225, and finally 300 mg.This case is also optimal, as can be seen from the graph of the cancer cells.The number of cancer cells decreased to zero in an approximation of a linear function relatively quickly.

Conclusions
In this study, we present a combination of two semi-analytical methods, HAM and MIM, to investigate the mathematical model of a BC treatment that combines an oral oestrogen receptor inhibitor called AZD9496 with palbociclib, a selective inhibitor of the cyclin-dependent kinases CDK4 and CDK6.A significant advantage of combining these two semi-analytical methods is that the mathematical model does not explicitly contain a small system parameter.By combining these methods, we revealed the hierarchy of the system of equations, which allowed us to divide the mathematical model into fast and slow subsystems.Thus, we studied only the fast system, while the slow system was frozen.In addition, we proposed two explicit analytical functions of the drugs, AZD9496 and palbociclib which are time-and dose-dependent functions of these drugs.This allowed us to control the dose and time intervals of treatment.Thus, we used different dose combinations of both drugs and found that there are two optimal combinations for the treatment of BC.One is that AZD9496 was not administered according to the standard protocol; that is, the dose of the drug changed cyclically.Additionally, the dose of palbociclib varied.The dose did not change regularly; that is, it increased and decreased.Second, the dose of AZD9496 constantly decreases while the dose of palbociclib constantly increases.

Figures and tables
In this section we present the values of the parameters used for numerical siulation apply to the mathematical model (1-11).In addition we present the solution profiles of the dynamical variables of the system.

Figure 1 .
Figure1.The solutions profiles of the cancer cells and the drugs protocols.The Palbociclib drug given is given arbitrarily i.e., the dose of the drug changes arbitrarily but the time intervals are constant.Whereas the AZD9496 drug is given according to the standard protocol, the dosage and time intervals are constant.

Figure 2 .
Figure2.The solutions profiles of the cancer cells and the drug protocols.The AZD9496 drug given is not according to the standard protocol, that is, the dose of the drug changes up and down, whereas the time intervals are constant.The dosage of the Palbociclib drug changes almost cyclic whereas the time intervals are constant.

Figure 3 .
Figure 3.The solutions profiles of the cancer cells and the drugs protocols.The dose of AZD9496 and Palbociclib drugs is constant in terms of the dosage and time intervals.

Figure 4 .Figure 5 .Figure 6 .
Figure 4.The solutions profiles of the cancer cells and the drugs protocols.The dose and the time intervals of Palbociclib are constant.Whereas the dosage of the AZD9496 drug changes.The dose does not change regularly, that is, it increases and decreases from the maximum 300 mg to the minimum of 50 mg.

Figure 7 .
Figure 7.The solutions profiles of the cancer cells and the drugs protocols.The dose of the drug Palbociclib is constantly decreased while the dose of the drug AZD9496 is constantly increased.Both of the drug's time intervals are constant.

Table 1 .
Parameter values using in the mathematical model (1-11).Day −1 pmol −1 Tumor growth rate induced by E2 Figure4presents the solution profiles of the cancer cells and drug protocols.The dose of palbociclib was constant in terms of dosage (125 mg) and time interval.