Analysis of a breast cancer mathematical model by a new method to find an optimal protocol for $HER2$-positive cancer

Abstract

Treatment of breast cancer (positive for HER2, i.e., ERBB2) is described by a mathematical model involving non-linear ordinary differential equations with a hidden hierarchy. To reveal the hierarchy of dynamical variables of the system being considered, we applied the singular perturbed vector field (SPVF) method, where a system of equations can be decomposed to fast and slow sub-systems with explicit small parameters. This new form of the model, which is called a singular perturbed system, enables us to apply a semi-analytical method called the method of directly defining inverse mapping (MDDiM), which is based on the homotopy analysis asymptotic method. We introduced the treatment protocol in explicit form, through an analytical function that describes the exact dose and intervals between treatments in a cyclical manner. In addition, a new algorithm for the optimal dosage that causes tumour shrinkage is presented in this study. Furthermore, we took the concept of protocol optimisation a step further and derived a differential equation that represents vaccination depending on tumour size and yields an optimal protocol of different doses at every time point. We introduced the treatment protocol in explicit form, through an analytical function that describes the exact dose and intervals between treatments in a cyclical manner. In addition, a new algorithm for finding the optimal dosage that causes tumour shrinkage is presented in this study. Additionally, we took the concept of protocol optimisation a step further and derived a differential equation that represents vaccination depending on tumour size and yields an optimal protocol of different doses at every time point.

Introduction

Breast cancer is a disease that involves a malignant tumour caused by uncontrolled growth of cells in a breast tissue.

Breast cancer is one of the most common cancers affecting women and is more prevalent after the age of 50. For instance, the American Cancer Society made the following estimates for breast cancer in the United States for 2018: $~$ 266,120 women will get a diagnosis of invasive breast cancer resulting in $~$ 40,920 deaths (American Cancer Society). There are many types of breast cancer such as ductal and lobular carcinoma in situ (non-invasive), invasive ductal carcinoma, and triple-negative breast cancer (invasive) (Breastcancer). The current study deals with $HER2$-positive invasive breast cancer (Krasniqi et al., 2019). Some females who get a diagnosis of breast cancer have a higher tumour level of $HER2$, which is an epidermal growth-promoting protein (epidermal growth factor receptor).

This kind of cancer tends to be more aggressive and to spread faster than other types of breast cancer (Palladini, 2009). Chemotherapy (Martine et al., 2005), hormonal therapy (Cancer.net and Hormonal Ther), and other therapies are used against $HER2$-positive breast cancer. In this paper, we discuss immunotherapy for $HER2$-positive breast cancer. Immunotherapy involves recruiting the human body's immune system to fight cancer cells. The immune system is capable of distinguishing self and non-self properties of cells and molecules and thus can eliminate pathogens and invaders. Although immunotherapy is based on the antigenicity of carcinogenic cells, and immune deficiencies are often accompanied by an increased risk of cancer, cancer cells are barely recognised as invaders because they are the body's own cells and their antigenicity is very low (Dunn et al., 2004). To overcome this challenge, scientists have developed a vaccine, called Triplex, that is expected to increase $HER2$ antigenicity and destroy the relevant cancer cells (Palladini, 2009), (De Giovanni et al., 2004). The Triplex vaccine includes three components: the $HER2$ gene product, which is the $p185$ oncoantigen (also known as $ERBB2$), $MHC$ class I molecules, and interleukin 12. Studies have shown that this type of vaccination via a specific chronic protocol (for life) can prevent mammary carcinogenesis and is helpful for the treatment (Palladini, 2009)- (Costa et al., 2017).

The incredible complexity of the immune system in general and its interaction with breast cancer in particular have been attracting many researchers who wish to understand the dynamics of this interaction. One of the main approaches to modelling this interaction is mathematical modelling; therefore, many fields have been engaged, such as biology, physiology, physics, and of course mathematics, for this task (Denise and Panetta, 1998). Mathematical models have been developed to describe the interactions among the immune system, breast cancer, and a relevant treatment and to investigate them for a better understanding of this disease and for designing better treatments (Botesteanu et al., 2016)- (Jorcyk et al., 2012). Hence, when one studies the phenomenon of cancer using mathematical models, with the aim of a better understanding of breast cancer progression and its antagonism with immunity, it is necessary to first model (in a reliable way) the biological phenomenon including translation of a complicated biological concept (such as death and interactions among the tumour, vaccination, and the immune system) into mathematical terms via a set of linear or non-linear partial differential equations (PDEs) or ordinary differential equations (ODEs) and second, to investigate this set of equations by analytical/numerical simulations and/or asymptotic methods. Subsequently, when we solve such systems and make a comparison with data obtained by $invivo$ or $invitro$ experiments, we have an opportunity to elucidate this phenomenon and the early stage of cancer progression and to predict its behaviour toward the vaccine and other immune-system entities, and we can formulate a putative optimal protocol of treatment.

Due to the complexity of solving the above type of models analytically, some $insilico$ models or simulators have been developed to interrogate and clarify the immune-system dynamics. For example, IMMSIM (National Institute for Cancer Research) is one of such computer programs for simulating the basic features and properties of the immune system on the basis of ODEs or cellular automata. One can use this software to model the response of a mammalian immune system to the presence of antigens and to predict its behaviour (Lollini et al., 2006). Another program is the SimTriplex simulator. This software also models the complex network of relations between the immune-system entities, and moreover, has been adapted specifically to mammary carcinoma and its antagonism with the Triplex vaccine (Pappalardo et al., 2005); accordingly, studies (Pappalardo et al., 2012)- (Bianca and Pennisi, 2012) on the modelling of this vaccine may include this simulator for a better explanation of $invivo$ experiments. Nonetheless, these computational methods do not allow for a qualitative, analytical, or asymptotic analysis. Another challenge with respect to the existing models (Botesteanu et al., 2016)- (Jorcyk et al., 2012) is stability examination. It is very difficult and sometimes even impossible to determine the stability of a model; this is a crucial step in researching an ODE system, especially when it represents a biological phenomenon and can have fateful implications for the patients.

In general, one can numerically solve an ODE system, even a partial one, of any order, for example, by means of different packages of Matlab, Mathematica 8.0, and/or Maple. Numerical solutions have drawbacks, but we will not address them in this study. One of the prominent disadvantages is that one can overlook a solution that goes off to infinity (for example, a singular point of the model). We propose to solve a given model semi-analytically/asymptotically, so that sometimes, we can obtain an analytical solution. The methods tested in this paper diminish the number of equations, and thus the reduced model can be investigated without losing information from the original model. Another significant advantage of the reduced model is run-time minimisation, which is, of course, an important feature.

To reduce the model under study, we apply a method called a singular perturbed vector field ($SPVF$). In a model described by a set of ordinary/partial, linear/non-linear differential equations, the rate of change of every dynamical variable in the system is hidden. The $SPVF$ method, in fact, enables us to expose this latent hierarchy. After implementing this method, we rewrite the model in new coordinates that practically represent a new version of a linear combination of the previous coordinates. Thus, we achieve decomposition of the system into two sub-systems: slow and fast. Hence, this approach allows for applying several asymptotic methods, for instance, perturbation methods that find an approximate solution to a problem (Winitzki, 2006). The homotopy analysis method ($HAM$) presented in ref (Liao, 2003). is an asymptotic semi-analytical technique for investigating mathematical models that are related to physics, biology, and other fields. Another useful asymptotic technique is the method of integral invariant manifold, which is a geometric asymptotic approach (Strygin and Sobolev, 1988). The main idea behind this method is to reduce the dimensionality of a system being considered, by splitting the system into fast and slow sub-systems and by studying these sub-systems on fast invariant manifolds. The above methods can be applied only to models that have the form of a singular perturbed system, i.e., the hierarchy of the model is predetermined. In particular, to use these methods, a given model has to explicitly include a small parameter(s), namely:

$$\begin{array}{l}{\mathrm{\epsilon }}_1{\dot{\overrightarrow{X}}}_1=F_1(x_1,\mathrm{...},x_k)\\ {\mathrm{\epsilon }}_2{\dot{\overrightarrow{X}}}_2=F_2(x_1,\mathrm{...},x_k)\\ .\\ .\\ .\\ {\mathrm{\epsilon }}_k{\dot{\overrightarrow{X}}}_k=F_k(x_1,\mathrm{...},x_k),\end{array}$$

where $\forall \{1<i<1\}$, ${\epsilon }_i\ll 1$ and ${\epsilon }_1<{\epsilon }_2<\dots {\epsilon }_i\ll {\epsilon }_{i+1}<\dots <{\epsilon }_k$. That is, there is a set of slowly changing variables in comparison to the other variables that change faster, and what divides the set of these equations is the fraction $\frac{{\epsilon }_i}{{\epsilon }_{i+1}}$.

In this paper, after exposing the hierarchy of the breast cancer model by using $SPVF$, we apply the method of directly defining inverse mapping ($MDDiM$) which is a semi-analytical technique presented in ref. (Liao and Zhao, 2016). The MDDiM is based on the $HAM$ method but without calculation of an inverse operator, as is typical in $HAM$; the latter calculation takes a lot of computing time. Another advantage of this method is that it is not necessary for the model to be presented in singular-perturbed-system form, i.e., the model can be presented without an explicit hierarchy.

The models that we investigated in this study are described in papers (Liao et al., 2014) and (Bianca et al., 2012). The model presented in ref (Liao et al., 2014). is a system of PDEs, and the model presented in ref (Bianca et al., 2012). is a system of ODEs. We combined these two models because we wanted to employ additional factors that influence the immune system and consequently affect the tumour. For the combined model, we proposed new functions for the vaccine dosage, two functions for the ODE model, and one function for the PDE model. Each function is defined in explicit form and describes the amount of a given vaccine for different combinations of protocols. We relied on the vaccination schedule that is composed of 4-week cycles of vaccination twice a week for 2 weeks, followed by 2 weeks of rest, and so forth (Bianca and Pennisi, 2012). Furthermore, the vaccination protocol is chronic, namely, the vaccination is administered throughout a patient's lifetime; accordingly, researchers in ref (Pappalardo et al., 2005). tried to develop a shorter protocol of vaccination without compromising treatment quality. We suggest an optimal protocol of vaccination based on different solutions of the $VC$ equation, yielding an optimal protocol. To the combined model, we add a new differential equation involving dynamical variable $VC$, i.e., variation of the vaccine amount with time t, depending on tumour size. The solution of this differential equation allows us to check the amount of the vaccine at any moment as a function of tumour size and thereby to control the vaccine dose.

To solve the combined model, we first solved the ODE model, obtained the dynamical variables of this model, and then substituted these variables into the PDE model. The purpose of this study is to first consider variables that are more dynamic among those related to the immune system and affect tumour growth. In addition, we propose a new protocol that varies depending on tumour size and hence is personalized.