Modeling of Mouse Experiments Suggests that Optimal Anti-Hormonal Treatment for Breast Cancer is Diet-Dependent

Estrogen receptor positive breast cancer is frequently treated with anti-hormonal treatment such as aromatase inhibitors (AI). Interestingly, a high body mass index has been shown to have a negative impact on AI efficacy, most likely due to disturbances in steroid metabolism and adipokine production. Here, we propose a mathematical model based on a system of ordinary differential equations to investigate the effect of high-fat diet on tumor growth. We inform the model with data from mouse experiments, where the animals are fed with high-fat or control (normal) diet. By incorporating AI treatment with drug resistance into the model and by solving optimal control problems we found differential responses for control and high-fat diet. To the best of our knowledge, this is the first attempt to model optimal anti-hormonal treatment for breast cancer in the presence of drug resistance. Our results underline the importance of considering high-fat diet and obesity as factors influencing clinical outcomes during anti-hormonal therapies in breast cancer patients. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-023-01253-1.


Introduction
Lifestyle factors such as age at menarche and menopause, body mass index, child birth and breast feeding, as well as genetic disposition, among others, are well-established breast cancer risk factors [1,2].However, much less is known about the role lifestyle factors play on breast cancer treatment response.Anti-hormonal treatment for estrogen receptor (ER) positive breast cancer constitutes a puzzling case in obese patients that requires more quantitative investigation.Approximately 75% of all breast tumors express ER, and most women with these tumors will receive anti-hormonal therapy [3].ER in breast cancer cells is activated by estrogen and it promotes cell proliferation and tumor growth [4].Anti-hormonal treatment with Aromatase Inhibitors (AI) decreases estrogen levels while anti-estrogen's block directly the action of steroids at the estrogen receptor [5].Interestingly, high Body Mass Index (BMI) and adiposity have a negative impact on AI efficacy [6][7][8][9][10][11][12][13][14].While the puzzle of the optimal anti-hormonal therapy in postmenopausal obese women is still unfinished, good monitoring of the suppression of estrogen levels with valid methods may guide treatment decisions during treatment with aromatase inhibitors [15,16].
An additional layer of complexity arises from the fact that ER-positive breast cancer cells may be resistant to anti-hormonal treatments.Resistance can arise due to multiple mechanisms that are not completely understood [17,18].Tumor cells can adapt to AI therapy after exposure for certain time (adaptive resistance), for instance due to the upregulation of ER expression or activation of alternative pathways conferring the cells survival and proliferative capacity.Instead, de novo or pre-existing resistance refers to the presence of estrogen independent cells before therapy.For instance, cells carrying specific mutations of the ER that confer constitutive ligand-idependent activity [19], which might lead to clonal selection under anti-hormonal treatment.The current paradigm consist of administering high AI doses to both obese and non obese patients, but this may not be the best strategy to avoid or delay drug resistance.
The aforementioned issues are difficult to quantify in preclinical and clinical settings and can benefit from more formal approaches.Here, we propose a new mathematical model, based on a system of ordinary differential equations (ODEs), to model the concentration of estrogen in the cancer tissue, which takes into account the local interplay between the tumor and fat tissues.We inform the model with data from mouse experiments that investigate the effect of obesity in breast cancer using two groups of mice, fed with control diet (CD) or with high-fat diet (HFD).Then, we incorporate AI therapy into the calibrated models, including de novo and adaptive resistance.To determine optimal therapeutic interventions in the CD and HFD cases, we formulate an optimal control problem (OCP) with the goal of minimizing the total tumor volume and the total amount of treatment that is used.We also compare the obtained optimal schedules with constant and alternating treatments.
Mathematical modeling of breast cancer dynamics under treatment have gained interest for long time [20][21][22][23][24][25][26][27][28][29].However, modeling of AI treatment in ERpositive breast cancer has received less attention so far.For example, an ODE model was proposed to understand pathway dynamics of ER-positive MCF-7 breast cancer cells under combination of Cdk4/6 inhibition and anti-hormonal therapies, including AI treatment [30].Similarly, Chen et al. proposed a mathematical model based on a system of ODEs to understand resistance to AI treatment driven by a shift from estrogen to growth factor receptors [31].In another study that uses stochastic differential equations and statistical physics techniques, the transitions under AI treatment between three different estrogen sensitive phenotypes were considered [32].To explain the dual effect of estrogen inducing both growth and regression of hormone-dependent breast cancer (referred as estrogen paradox), Ouifki and Oke proposed an ODE model and determined conditions to eliminate cancer recurrence for long-term treatment based on stability analysis [33].Cancer treatment scheduling optimization by means of OCPs has received considerable attention [34][35][36][37].For instance, OCPs were proposed to optimise treatment schedules of chemotherapies [38][39][40], angiogenic inhibitors [41], cytotoxic and antiangiogenic therapies [42], immunotherapy via a dendritic cell vaccine [43] and combination therapies [44,45].In addition, resistance to chemotherapy [46,47] or combination of chemotherapy with ketogenic diet [48] were also studied using OCPs.Another recent study investigated the optimal combination of doxorubicin and HER2 targeting drug trastuzumab, for a murine model of human HER2 positive breast cancer [49].To the best of our knowledge, the present work is the first modeling study to account for anti-hormonal treatment using AIs in the presence of drug resistance in an optimal control framework.
The paper is organized as follows: In the following section 2, we formulate the dynamical model, prove some basic properties and we proceed with model calibration.Sec. 3 is dedicated to the model extension for AI treatment and resistance.In Sec. 4, we formulate the OCP and derive the optimality system.
Then, we proceed with results in Sec. 5 to compare various relevant scenarios for anti-hormonal treatment.We conclude by discussing the main conclusions, limitations and future directions.

Mathematical model development and calibration
In this section, we develop a basic ODE model for the interaction of ERpositive breast tumor cells, estrogen and fat for the postmenopausal situation, and demonstrate some useful mathematical properties of its solution.Our model can describe the contribution of fat intake differences to estrogen and tumor growth over time.We have been inspired by the mice experiment conducted by Hillers et al. [50] comparing tumor growth in CD and HFD mice, and we use the data obtained in that study to bring our model closer to reality.We proceed by describing the experiments and available data that inspired our model.We then present model equations and the assumptions they are based on.Then, we discuss mathematical properties of the model.Lastly, we explain the details of the model calibration.

Experimental data
Hillers et al. investigated the influence of obesity on breast tumor size and stromal cells within the mammary adipose tissue [50].We use data from that study that utilized breast cancer cell line EO771 derived from a spontaneous mammary adenocarcinoma from a C57Bl/6 mouse.EO771 cells are considered to be a model of luminal B breast cancer subtype and are known to respond to anti-estrogens [51].Specifically, mice were fed with CD (10% kcal from fat, Test Diet 58Y1) or HFD diet (60% kcal from fat, Test Diet 58Y2).A total of 1 × 10 6 EO771 tumor cells were mixed with 2.5 × 10 5 adipocytes taken from CD or HFD mice.After pelleting this mixture of cancer and fat cells, it was injected bilaterally into the inguinal mammary glands of 8-week-old female mice fed with CD.In total, we have the data of eleven mice, five where the fat cells come from mice fed with CD and six where fat cells come from mice fed with HFD.For each tumor independently, tumor volumes were measured at days 10, 13 and 15, as depicted in the study [50, Fig. 2B].In addition, the number of adipocytes at day 15 was quantified, see [50, Fig. S2F], and we use it to estimate the amount of fat in the tumor tissue.Details are provided in the Supplementary Material 7. Mice were euthanized when tumor reached the humane endpoint of 15 mm in diameter.In summary, we obtained the initial conditions and six independent measures of tumor volume at three time points (days 10, 13 and 15) for each condition (CD and HFD), and two data points for fat volume at day 15, one for each condition (CD and HFD).

Model development
In order to quantify the effect that the fat-induced production of estrogen has on tumor growth, we model the temporal dynamics of tumor volume T := T (t) (mm 3 ), estrogen concentration E := E(t) (pg/g) and fat volume F := F (t) (mm 3 ) in the tumor tissue at time t (days).The model is based on the following six assumptions: 1. Tumor volume follows logistic growth [52] .
2. Tumor growth rate depends on the estrogen level [51].
3. Fatty tissue is the major source of estrogen in the tumor [53].Circulating estrogen concentrations are proportional to adipose mass in postmenopausal women [54], so we assume that estrogen is produced by fat at a constant rate.4. Estrogen is washed out from the tumor tissue at a constant rate [55]. 5. Tumor cells use fat as an energy resource [56,57].6.While mice are fed CD for 15 days, there is no growth source for fat volume.
Therefore, diet-based difference in fat volume (CD vs HFD) are accounted for simply due to the amount of fat volume at day 0. 7.In the experiments that we model, most implanted adipocytes survive.Cancer cells are known to produce growth factors and cytokines that support the survival of adipocytes.Therefore, we did not include fat decrease due to adipocyte death.
A flow diagram depicting the interactions between model variables T , E and F is presented in Fig. 1(a).
Consequently, we propose the following system of ODEs: dT dt change in tumor volume

=
prolif eration rate triggered by estrogen The parameters k 1 , a 1 , m 1 , r, µ, α and initial conditions T 0 , E 0 and F 0 are all non-negative real numbers.Eq. (2.1a) represents the tumor logistic growth, where the growth rate is assumed to follow Michaelis-Menten kinetics, g(E) = k1E a1+E .Parameter k 1 is the maximum growth rate for high estrogen levels and a 1 is the estrogen concentration at which the growth rate is half-maximum.Estrogen is produced by fat but this production is inhibited by AIs.Estrogen is also naturally washed out.In both diagrams, the lines ending with an arrow represent positive feedback whereas the lines ending with a bar denotes negative feedback.
Parameter m 1 is the inverse carrying capacity of the tumor.Eq. (2.1b) models the change in estrogen concentration when it is produced by fat at a rate r and washed out from the tumor tissue at a rate µ.The last equation (2.1c) accounts for fat consumption by tumor cells at a rate α.The values of these parameters are not known and will be estimated from data.

Model properties
Next we prove that the solution to model (2.1) exists, it is unique, non-negative and bounded.These properties will be used later.
Proposition 1 Eq.(2.1) with non-negative initial conditions has a unique solution that is non-negative and bounded from above for all t ≥ 0.
Proof As the right-hand side of the model (2.1) and their partial derivatives are continuous on R×R 3 , it follows from the Cauchy-Lipschitz theorem that the existence and uniqueness of the solution are guaranteed [58,Ch.15].
To prove that the solution to Eq. (2.1) is non-negative for all t ≥ 0, we use the method of separation of variables.Firstly, Eq. (2.1c) leads to Since F and r are non-negative, we can rewrite Eq. (2.1b) as Thus, E ≥ 0. Finally, Eq. (2.1a) leads to Therefore, T ≥ 0, E ≥ 0 and F ≥ 0 for all t ≥ 0.
To prove that the solution to Eq. (2.1) is bounded from above, we observe from Eq. (2.1a) that Thus, the solution is bounded from above.□

Model calibration
In this section, we make use of the experimental data described in Sec.2.1 to inform our basic model.Direct measures of model parameters are not available in this experimental setup, and we do not have enough data to do formal statistical inference for all parameters.The main argument we used to fix some parameters was identifiability of the remaining free parameters.To that goal, we decided to fix three of them to reasonable values and made extra assumptions to fix the initial estrogen concentration and fat volume.We then use the available data to calibrate the rest of the parameters for which we lack any information and show that the problem is practically identifiable by using profile likelihood [62].
The initial amount of fat in the tumor tissue was measured only at day 15.For simplicity, we assume that the level of fat under CD stays constant and it has not changed since the beginning of the experiment (see [50,Fig. S2F]).We acknowledge this is a limitation and an initial fat measurement would have made our results more solid.As estrogen is mainly produced by fat, we also assume that estrogen concentration is proportional to fat volume at baseline.Indeed, estrogen concentration in mice under different, but comparable, conditions was measured between 150-1500 pg/g [59, Fig. 2].For estrogen concentration in our model to lie within those measures, we assume that the ratio of estrogen concentration to fat volume is around 3.4.
We then find reasonable values for parameters m 1 and µ.We obtain the half-life of estrogen in breast tumor tissue from [55], t 1/2 =2.8 h.Therefore, µ, that represents the washout rate of estrogen from tumor tissue, can be computed as µ = ln(2)/t 1/2 = 0.25 h −1 = 5.94 day −1 .In the experiments, mice were euthanized after the tumor reached 15 mm in diameter, corresponding to a volume of approximately 1767mm 3 assuming a spherical tumor.We simply set m 1 = 1/2000 mm −3 equally for CD and HFD as a larger value than the highest tumor volume in the data set.Fixing m 1 and µ still did not solve the nonidentifiability problem, but we discovered that fixing r in addition solved this issue.Parameter r is fixed as 20 pg/g mm −3 day −1 based on the assumption that estrogen is at the steady-state in the beginning of the experiment which leads to E = rF µ .Estrogen concentration roughly satisfies 150 ≤ E = rF µ ≤ 1500 [59, Fig. 2].By multiplying both sides of this inequality by µ = 5.94, we reach 891 ≤ rF ≤ 8910.We divide both sides by F CD (0) = 50 and F HF D (0) = 360 separately, that results in two inequalities 17.82 ≤ r ≤ 178.2 and 2.475 ≤ r ≤ 24.75.The intersection of these inequalities gives a range for the parameter r which is 17.82 ≤ r ≤ 24.75.Therefore, we simply choose r = 20.
We perform model calibration and profile likelihood calculations in Data2Dynamics [60,61].We fixed the lower and upper bounds for the parameters as 10 −7 and 10 4 in the optimization problem, respectively.Based on the method of profile likelihood [62], our model is practically identifiable (See.7).We list the obtained parameter values in Table 1.We perform sensitivity analysis for all the parameters in Sec. 7.
Fig. 2 shows the simulation results for CD (left panel) and HFD (right panel) obtained with the parameters in Table 1.The y left-axes correspond to tumor or fat volume, whereas the y right-axes denote the estrogen level.Data points with error bars for tumor and fat volume are marked with red circles and black crosses, respectively.We observe that simulation results for tumor and fat volume agree well with the data, whereas estrogen level stays within a biologically meaningful interval.We observe that the initial estrogen level and fat volume are higher for HFD than CD (see, Table 1.).Temporal evolution of estrogen level and fat volume is similar for both diet types, since fat is assumed as the source of estrogen production.Indeed, tumor associated with HFD increases faster than for the other, due to more estrogen release from HFD fat volume.Original model (2.1) expresses the changing dynamics of tumor volume, estrogen level and amount of fat in case of no treatment and simulation results agree well with the available data.The next step is to extend this model to account for AI treatment by considering sensitive and resistant tumor subpopulations.In this way we will be able to study drug resistance to endocrine therapy for ER-positive breast cancer.

Model extension for resistance to aromatase inhibitor treatment
Aromatase inhibitors, despite of being an effective treatment choice for ERpositive breast cancer, may suffer from drug resistance [18,63].To investigate different treatment schedules, including constant, alternating and optimal antihormonal treatment, we consider tumor heterogeneity in terms of sensitive and resistant subpopulations under the following assumptions: 1. Breast cancer cells are either sensitive or resistant to estrogen deprivation with AIs.In reality, there could be more than two tumor subpopulations, since development of resistance is considered as a progressive mechanism and cells may shift form one stage to another over time [64].For simplicity, we assume that there are only two tumor subpopulations, called sensitive and resistant.2. Both sensitive and resistant populations follow logistic growth [52].3. Growth of sensitive cells is triggered by estrogen [65].4. Sensitive cells die under low estrogen concentrations [65]. 5. Sensitive cells adapt to low estrogen levels and become resistant cells [31,32].6. Resistant cells do not die under low estrogen concentrations [31,32].7. Fat volume follows logistic growth [66].8.Both sensitive and resistant cells consume fat as energy resource [66].
Consequently, Eq. (2.1) is extended with the sensitive cell population S := S(t) (mm 3 ) and the resistant cell population R := R(t) (mm 3 ) in the tumor tissue at time t (days): ) with the non-negative initial conditions S 0 , R 0 , E 0 and F 0 .Eq. (3.1a) expresses the logistic growth of sensitive cells over time together with death and adaptation terms.Parameter m 1 is the inverse of maximum tumor size and η is the competition parameter scaling inhibition of sensitive cells' growth by resistant cells.Sensitive cells die if estrogen level is smaller than a 2 while they adapt to estrogen level below a 3 and become resistant.Parameter c is the maximum death rate and l denotes Hill's coefficient.Eq. (3.1b) models evolution of resistant cells with the growth rate k 3 .Eq. (3.1c) stands for dynamics of estrogen level where the parameter p, 0 < p ≤ 1, reduces the effect r to p • r due to aromatase inhibitors.Eq. (3.1d) models the change in fat volume with logistic growth so that effect of fat growth to anti-hormonal treatment could be investigated.Parameters k 2 and m 2 are the growth rate of fat and inverse carrying capacity of fat, respectively.The carrying capacity could model how the body is prone to accumulate fat depending on the life style or other metabolic conditions.In addition, both sensitive and resistant cells consume fat as energy resource at the rate α.A diagram depicting the interactions between the extended model variables S, R, E and F is presented in Fig. 1(b).
We assume that the parameters that we calibrated in the basic model are not affected by the treatment and we use them in the extended model.As we do not have data under treatment, we explore the effect that new parameters have by testing different values.

Model properties
Proposition 2 Eq.(3.1) with non-negative initial conditions has a unique solution that is non-negative and bounded from above for all t ≥ 0.
Proof Existence and uniqueness of the solution is standard and analogous to Proposition 1.Thus, we prove here that the solution to Eq. (3.1) is non-negative and bounded from above for t ≥ 0. Similar to Theorem 1, we can prove that E, F ≥ 0 for t ≥ 0 by the variation of constants formula.For Eq. (3.1a), we have Since S ≥ 0 for t ≥ 0, Eq. (3.1b) can be written as Then, we get We can prove that F and E are bounded from above similar to Theorem 1.On the other hand, using Eq.(3.1a)-(3.1b),we obtain the sum Since S and R are non-negative, it means that S and R are bounded above.Then, we complete the proof.□

Treatment modelling
We will investigate differences between constant, intermittent and optimal anti-hormonal treatment.Constant treatment is implemented through the parameter 0 ≤ p ≤ 1 in Eq. 3.1.The value p = 1 corresponds to no estrogen deprivation treatment and smaller values of p models AI treatment with inhibition of estrogen production.
Alternating treatment refers to a pre-scheduled treatment scenario with u I := u I (t) and it is implemented by modifying Eq. (3.1c) to where In the next section, an OCP is constructed to investigate the optimal value of p as a time-dependent function, and results obtained with the optimal endocrine therapy are compared with the constant and alternating treatment.

Optimal control problem for anti-hormonal treatment
We aim to investigate optimal AI treatment schedules that minimize the total number of cancer cells together with the pharmaceutical intervention over a prespecified time interval [t tr , t f ].We do not include an equation representing the drug as often done for optimal chemotherapy scheduling in the literature (see, for example, [67,68]).Instead, we model the effect of AI treatment through a continuous control function u := u(t).AIs act by lowering the estrogen production, so we replace the parameter p in Eq. (3.1c) by the function 1 − u.
We formulate the OCP as follows: minimize the cost functional subject to ) ) ) where, Our aim is to find an optimal control u * such that J (u * ) = min u∈U J (u).
We note that constructions of linear or quadratic cost functional in the control function u results in not only biologically but also mathematically different interpretations.While quadratic OCPs have a single extremum and result in continuous controls, linear OCPs result in bang-bang controls and mathematical analysis becomes more complicated due to singular or bang-bang controls that result in non-differentiable solutions curves.We refer readers to the following studies for a detailed comparison [69][70][71].In addition, the parameters ω S , ω R and ω u in Eq. (4.1) can be set to balance the size of the different terms.
In the present case, inclusion of the term u 2 in the Eq.(4.1) is justified by the treatment side effects.Side effects of AI include from hot flushes to cardiovascular events, vaginal bleeding and bone loss [72][73][74].Our quadratic choice reflects the fact that the increase in side effects is negligible for small amounts of therapy and that side effects increase as function of u, rather than increasing at a constant rate as in the linear control.
Theorem 1 There exists an optimal control u * with a corresponding solution (S * , R * , E * , F * ) to the model (4.2) with non-negative initial conditions that minimizes (4.1) over U.
Proof The proof is based on several steps according to the study of Fleming and Rishel [75,Corollary 4.1].Firstly, we observe that the coefficients in Eq. (4.2) and its solution are bounded on a finite time interval, so the admissible control set U and the corresponding solution with initial conditions are non-empty [76, Thm 9.2.1.].Secondly, the admissible control set U is closed and convex.In addition, the right-hand side of the system (4.2),namely ⃗ f (t, ⃗ X, u) with ⃗ X = (S, R, E, F ) T , is continuous, since the system has positive parameters and the non-negative solution by Proposition 2. Indeed, it is bounded above by a linear combination of the bounded control and the state as due to bounded solution (by Proposition 2) and positive parameters in the model for some positive constant C. On the third line, we use the relation The integrand of the objective functional is convex on U due to the quadratic term.Indeed, it is bounded as with some positive constant Ĉ.Thus, we can conclude that an optimal control u * exists.□ Theorem 2 Given an optimal control u * and solution to the system (4.2) for (4.1), there exist adjoint variables λ i := λ i (t) for 1 ≤ i ≤ 4 such that ) ) ) with Proof Following references [67,77], the Lagrangian is constructed as where the Hamiltonian H is defined as and ξ i (t) ≥ 0 are penalty multipliers such that From the Pontryagin's Maximum Principle, we can derive the adjoint equations by obtaining partial derivative of the model (3.1) with respect to S, R, E and F , respectively.Indeed, we get with To obtain an expression of the control, we differentiate the Hamiltonian with respect to u as and project it onto the admissible set of controls.□

Implementation of the optimal control problem
Third-generation AIs (anastrozole, letrozole and exemestane) reduce wholebody aromatisation by >90% (summarised in ref. [78]).However, limited local estrogen production in tissue compartments cannot be totally ruled out.Also, it is possible that some cells could locally produce some estrogen under treatment [79,80].Therefore, we assume that the maximum drug dose does not eliminate the total estrogen in the vicinity of tumor.This could be done simply by setting a threshold value on the control function.We use u a = 0 and u b = 0.99, where u a corresponds to the case of no treatment and u b refers to the strongest possible treatment.
The optimality system consisting of the state equation (4.2), the adjoint equation (4.7) and the optimality condition (4.8) form a nonlinear system of equations, so we obtain the numerical solution via forward-backward sweep (FBS) method [81].As explained by Lenhart and Workman [81], the FBS method requires initiation of a feasible control function to solve the state equation forward in time.Then, the adjoint equation is solved backward in time and the optimality condition is updated at each iteration until the stopping criterion is satisfied.Here, the idea is to find a feasible optimal control iteratively.The update strategy of the control could be done in different ways such as taking average of the current (u cur ) and previous control (u pre ) or their convex combination [81].Here, we apply the approach "greedy" convex combination studied by Vatcheva et al. [82,Sec. 3] to cover a large range of control combinations during optimization and avoid stagnation."greedy" convex combination refers to expressing the control as u s = (1 − s)u pre + su cur where s ∈ (0, 1) is selected in such a way that the smallest value of J (u s ) is achieved in that iteration.The parameter s is not fixed as opposed to the averaging or convex combination, it may vary in each iteration instead.The stopping criterion in this paper is based on the relative error of the current and previous state, adjoint and control functions.The program is terminated when a relative error less than 10 −5 is achieved.
Simulations in this study were performed using MATLAB ® R2022 [83].We used ode15s solver to obtain the numerical solution of the differential equations and fmincon function in the model calibration step.All data and code are available (see data and code availability part for the details).

Simulation results
We focus in simulations of the extended model (3.1) that explore the effect of threshold values a 2 and a 3 .These two values correspond to the estrogen concentrations below which cancer cells die or become resistant, respectively.Thus, simulation scenarios using different threshold values represent treatment in hypothetical tumors with differential sensitivities and rates of resistance to the local estrogen availability.For each case, we simulate three different treatment types: constant treatment, alternating treatment and optimal antihormonal treatment.We use the parameter values which are common in both the first and the extended model.For the others, we either fix their values or explore their impact in simulations.We list all parameter values in Table 2.
For constant treatment, the parameter p is chosen from the set {1, 0.025, 0.0125, 0.01, 0.001}.For alternating treatment, we set u b = 0.99.Treatment is started on the date corresponding to the earliest time point t := t tr at which S + ηR < 1 4m1 so that the tumor reaches a detectable size.The final simulation time is fixed as t f = 25 to obtain a unique optimal control (see [84,Sec. 4], for detailed discussion).
We choose the weight coefficients ω S , ω R , ω u in the cost functional as one or hundred to model different penalization strategies.For instance, values ω S , ω R > ω u refers to penalization of tumor cells more than treatment cost. 1 See appendix for computational details.
5.1 Scenario I with a 2 = 20 pg/g, a 3 = 1 pg/g, k 3 = k 1 2 5.1.1Scenario Ia: Adaptive resistance with R 0 = 0 We first explore a scenario without preexisting resistance, where initially all cells are assumed to be sensitive to treatment and no resistant cells exist.Instead, endocrine resistance may arise due to adaptation to low estrogen levels.Tumor cells may die if estrogen level is below a 2 = 20 pg/g and they may become resistant if estrogen concentration falls below a 3 = 1 pg/g.Proliferation rate of the sensitive cells is assumed to be equal to half of the growth rate of sensitive cells.In Fig. 3, we show response to constant treatment by displaying the change in tumor size for different values of p.The solid line corresponds to the case where no treatment is applied, i.e., p = 1 and tumor reaches to the carrying capacity as time passes.We mark in the figures the time point at which treatment is started with a dashed vertical line and we observe that treatment is started earlier for HFD than CD.We observe that anti-hormonal treatment results in eradication of tumor for the values p = 0.0125 and 0.01 for both CD and HFD, whereas the case p = 0.025 does not lead to tumor eradication for HFD.In case of a drug inhibiting estrogen production 99.9%, i.e., p = 0.001, drug resistance is observed.We mark the time point at which treatment is started with a dashed vertical line.
Next, we proceed with alternating treatment.We simulate two different schedules with long and short treatment phases in Fig 4 .The solid curve refers to the tumor size without treatment.For CD, alternating treatment with shorter phases causes tumor volume to stay within a range.Instead for HFD, tumor size grows over time with respect to the baseline tumor volume.On the other hand, treatment with longer phases leads to tumor reduction for CD, while it stays within a range for HFD.
In case of optimal treatment scheduling, we observe in Fig. 5 that the tumors are eradicated for both CD and HFD (solid line for no treatment, dashed line for optimal treatment with ω R = ω S = ω U = 1).For comparison, optimal control functions u(t) are shown in Fig. 6 for different values of the weight constants.The case ω R , ω S > ω U leads in both CD and HFD to maximum treatment for almost the entire studied period, whereas treatment could be stopped earlier if ω R , ω S < ω U .In other words, penalizing tumor cells more than treatment results in longer treatment.There is no big difference between control functions in terms of diet, except for a slightly larger duration of treatment for HFD.To observe the effect of treatment in more detail, we present in Fig. 7 the dynamics of all variables in the case ω R = ω S = ω U = 1.In this example, optimal treatment maintains the estrogen level between a 3 and a 2 , so sensitive cells die, but no resistance occurs.This happens in spite of an increasing fat volume.Therefore, optimal treatment results in successful elimination of tumor without causing drug resistance.

Control function
Fig. 6: Scenario Ia: Optimal control function u over time t for three different combinations of weight coefficients ω R , ω S , ω U .Dashed and dash-dotted curves refer to the optimal treatment schedules for CD and HFD, respectively.Solid line denotes the maximum treatment.We mark the time point at which treatment is started with dashed and dash-dotted vertical lines for CD and HFD, respectively.Next we investigate the influence of a preexisting resistant sub-population on the success of constant, alternating and optimal anti-hormonal treatment schedules.In this case, endocrine resistance arise by clonal selection of cells that are endocrine independent for some reasons.In Fig. 8 we show response to constant treatment by displaying the change in tumor size for different values of p.We observe that constant treatment is unsuccessful to eliminate the tumor.We plot the change of tumor size over time for short and long alternating treatment phases in Fig. 9 where 75% of the cells are sensitive and 25% are resistant at the beginning of the simulation.For shorter drug holidays, tumor size increases compared to the initial tumor size.For longer on-off periods, the tumor volume remains within a bounded range for both CD and HFD but the oscillations between remission and growth are bigger in the HFD case.Moreover, the final tumor volume is larger for both cases in comparison with the case of no preexisting resistance showed in Fig. 4. We mark the time point at which treatment is started with a dashed vertical line.
We plot the results obtained with optimal treatment in Fig. 10.While treatment decreases the tumor volume in both CD and HFD cases, resistance cells proliferate and drug resistance occurs.We present the temporal evolution of the model variables in detail in Fig. 11.We can see that sensitive cells are killed but resistant cells increase in size as a result of drug-resistance.Optimal control profiles are similar to the case of adaptive resistance in Fig. 6, so we do not present it here.5.2 Scenario II: Adaptive resistance with a 2 = 10 pg/g, Here we investigated a situation with no preexisting resistant cells but where the estrogen thresholds for which the cells die or are converted to resistant are closer than in the previous case.Fig. 12 shows the case for constant treatment.All tested treatment cases, except p = 0.001, result in tumor elimination for CD, but p = 0.01 leads to tumor reduction until day 25 and higher values of p suppresses tumor growth for HFD.Alternating treatment instead, leads to oscillations in tumor size but with a decreasing trend for CD, whereas for HFD a sharp increase in tumor population is observed during drug holidays (see Fig. 13).
The results for optimal treatment shown in Fig. 14 reveal that tumor is eradicated for CD, similar to Fig. 5. Interestingly, the tumor remains at the end of the treatment for HFD.We compare optimal treatment schedules in Fig. 15 for CD and HFD.We note that treatment must be applied for long time for HFD than CD, while it could be relaxed earlier for CD.Thus, optimal anti-hormonal treatment gives the most promising results among three different treatment choices in terms of reduction in tumor volume and time to lessen treatment could also be seen.Finally, we investigate a scenario where death and conversion terms are equivalent, namely a 2 = a 3 = 10.We present temporal evolution of all model variables for constant treatment in Fig. 16.Estrogen level is successfully decreased, but it leads to drug resistance for CD for all choices of the parameter p.On the other hand, for HFD, the case p = 0.025 is not strong enough to kill sensitive cells, so resistance cells do not proliferate.However, other treatment choices result in resistance and treatment fails.On the other hand, alternating treatment is not a successful strategy (see Fig. 17).We mark the time point at which treatment is started with dashed and dashdotted vertical lines for CD and HFD, respectively.
Finally, optimal AI treatment results in drug resistance as seen in Fig. 18.Initial tumor size reduction is followed by cell proliferation.Even though treatment is stopped earlier for CD than HFD (see Fig. 19), it is not possible to eliminate resistance due to equal cell death and conversion terms in the model.
A detailed picture of model variables is presented in Fig. 20 and it reveals that treatment kills sensitive cells due to low estrogen level; but, then resistance occurs.

Conclusions of the simulation results
We compared outcomes for different treatments in a series of hypothetical tumors with differential sensitivities and rates of resistance to the local estrogen availability.We observed that in tumors where the difference between estrogen thresholds for cancer cells to die and to adapt to low estrogen levels is large, then constant treatment with an appropriate dose or optimal treatment Optimal anti-hormonal treatment for breast cancer   are the best for the case of only adaptive resistance.However, if the difference between the thresholds is smaller, then optimal treatments are better, specially in the HFD case.In case of preexisting resistance, thresholds for cancer cells to die and to adapt to low estrogen levels is large, optimal treatment or constant treatment with appropriate dose gives the best outcome.When death of cancer cells and their adaptation to level of estrogen occur at the same threshold value, optimal treatment is best choice.Importantly, treatment outcome and optimal treatments schedules differ based on diet.We mark the time point at which treatment is started with dashed and dashdotted vertical lines for CD and HFD, respectively.

Discussion
Given the rising obesity rates around the world, novel strategies are urgently needed to evaluate and optimise endocrine treatment of breast cancer in women with high BMI.In this study we focused on modeling the effect that fatinduced production of estrogen has on tumor growth.While our model is able to capture the trends in the experimental data for CD and HFD mice, we recognise that other factors associated with the adipose tissue and not considered in our current model, such as inflammatory cytokines, leptin or insulin, could be influencing tumour growth differently in the CD and HFD cases.These are subjects that deserve further investigation [85].By incorporating AI treatment and resistance in our model, we can simulate treatment outcomes in CD and HFD mice.However, as we do not have data on treatment, the choice of parameters related to sensitivity and resistant to treatment were made by explorative simulations.For instance, we assumed cost of resistance in the sense that the growth rate of resistance cells is smaller than the growth rate of sensitive cells.Otherwise, rapidly increasing resistant cells would always dominate the tumor.In addition to this, more than two tumor subpopulations with differential drug-response to AI could exists.When AI treatment data in these mice are available, it would be possible to obtain the number of subpopulations, their fractions and their growth rates trough a novel phenotypic deconvolution method [86].
Besides constant and alternating treatments, we investigated optimal scheduling trough OCPs.In this framework, we underline that one of the theoretical challenges is to prove uniqueness of the optimal control on a specific time interval [0, t f ], since the value of t f cannot be found explicitly, and it is bounded by some constants depending on the solutions of the state and adjoint equation.We observe that a larger time interval leads to convergence issues, which is an indication of the uniqueness of the solution on a smaller time interval.We have also experienced that the more complicated the ODE model used in the OCP constraint is, the smaller the time interval where a unique solution can be found.Furthermore, uniqueness could be proved using constant tumor growth rate, but we believe this is not a correct representation of ER-positive tumor subtype.
Being a breast cancer modeling study with optimal control analysis, Oke et al. constructed a model of four variables (including normal cells, tumor cells, natural killer cells and estrogen concentration) with implementation of anticancer drugs and a ketogenic diet [48].They modelled the ketogenic diet as a parameter affecting tumor growth, while anti-cancer drug was modeled as an intervention strategy leading to tumor death, and estrogen concentration to decrease, so that suppression of immune cell activation was relaxed.In addition, optimal values of the parameters corresponding to anti-cancer drugs and ketogenic diet were searched to minimize the total tumor size and estrogen concentration on a prespecified time interval within a quadratic optimal control setting.The authors noted that activities of cancer cells are reduced with the introduction of a ketogenic diet and they underlined the risk of ketoacidosis as a results of too much ketogenic diet.The authors found based on stability analysis of tumor-free equilibrium point that tumor cells could be eliminated with treatment and ketogenic diet, if the reproduction number of the system was reduced to a value less than one.This is in contrast with our simulations, where HFD does not result in better treatment outcomes.Interestingly, it has been shown that different fat diets, i.e. based on olive vs corn oil, influence breast tumor growth and progression differently [87,88], adding complexity to the challenge of optimizing breast cancer treatment and diet.
Overall, the most striking observations from our simulations are that optimal aromatase inhibitor treatment schedules and the corresponding outcomes differ based on diet, which suggests that low fat diet and other measures to reduce the amount of fat could be introduced to improve treatment outcomes in obese patients.In our ongoing studies, we are modeling such patient-specific treatments making use of individual level data from the NeoLetExe trial [89], a neoadjuvant study exploring the lack of cross-resistance between the aromatase inhibitor letrozole and the aromatase inactivator exemestane.The effect of switching to a low-fat diet is not necessarily immediate, because it also depends on the lifestyle and how the body is prone to accumulate fat.Although our extended model might be able to capture lifestyle effects different than diet, this remains to be investigated.Optimal anti-hormonal treatment for breast cancer Data and code availability Simulations in this study were performed using MATLAB ® R2022 [83].All data and code used in this article are publicly available in the online repository of the Oslo Center for Biostatistics and Epidemiology (OCBE).

Supplementary Material
Computation of fat volume at t = 15 and carrying capacity of adipocytes: Growth of adipose tissue is based on two processes affected by genetic and diet differences, namely cell number increase (hyperplasia) and cell size increase (hypertrophy) [90].Therefore, computation of the amount of fat at days t = 0 and t = 15 in the system are based on these two mechanisms.We start by calculating the average number of adipocytes per image and per mm 2 for CD and HFD, separately, based on the study [50, Fig. S2F].
On the other hand, Bozec and Hannemann report that adipocyte size increases in mice fed with HFD, so we set the diameters of adipocytes for CD and HFD as 0.1 mm and 0.2 mm, respectively, as reported in the study [91,Fig.6C],which leads to two different values for the amount of fat at time t = 15 for CD and HFD.In addition, mice are fed with control diet after tumor injection in the experiment, so we assume that carrying capacity of adipocytes associated with mice fed with CD and HFD are determined by the amount of fat for HFD at t = 15.
Procedure to estimate the amount of fat in the cube with volume V could be summarized as follows: Let n be the number of adipocytes per mm 2 and d be the average diameter (mm) of an adipocyte.Then, 3  √ V is the size of the cube in mm with volume V and V ×n adipocytes per mm 2 .Finally, the approximate number of adipocytes in the cube with volume V could be aproximated via the product .
Profile likelihood calculations We present profile likelihood calculations for model (2.1) after performing calibration [62] in Fig. 21.We used "arPLEInit" and "ple" functions of d2d software.The figure indicates practical identifiability of the model corresponding to the 95% confidence level for the parameters a 1 and k 1 .For parameter α, this threshold is above 83%.

Sensitivity analysis
We present the results of the global sensitivity analysis for Eq.2.1 to capture the relative changes of the variables with respect to the parameters at days 5, 15 and 25.Here, the method of partial rank correlation coefficient (PRCC) presented by Marino et al. is used [92].
Global sensitivity analysis requires a uniformly distributed sample space which is constructed for each parameter with 1000 sample values.Here, intervals for each parameter are constructed with the end points as ± 5 % of the baseline parameter.LHS-PRCC MATLAB code given in the web site [93] have been used and modified for the current model.The idea behind PRCC is to assign a value between -1 and +1 to each relation.The magnitude of this value determines the strength and the sign indicates the trend of the relation between the parameter and the variable.We briefly compare the most sensitive parameters where sensitivities for CD and HFD are grouped on the left and right panel in Fig. 22.Firstly, we observe that parameters k 1 and r have positive sensitivities for tumor volume, whereas negative sensitivity with respect to µ decreases over time in magnitude.Sensitivity of estrogen concentration is positive for r and negative for µ, as expected.In addition, fat volume is negatively sensitive to k 1 and positive sensitivity with respect to m 1 increases over time and it is effect is smaller for HFD than CD case.

Fig. 1 :
Fig. 1: (a) Modelled interactions between the volume of tumor cells, T , fat volume, F , and estrogen concentration, E, in Eq. (2.1).Tumor volume grows as a consequence of cancer cell proliferation, which is triggered by estrogen.Estrogen is produced by fat and it is washed out.Tumor cells consume fat as energy resource.(b) Modelled interactions between the volume of resistant cells, R, volume of sensitive cells, S, fat volume, F , and estrogen concentration, E, in Eq. (3.1).The volume of both sensitive and resistant cells grows as a consequence of cell proliferation.While the growth of sensitive cells is triggered by estrogen, that of resistance cells is estrogen independent.Sensitive cells can die under the influence of estrogen or can adapt to low estrogen levels and become resistant.Both sensitive and resistant cells consume fat.Fat volume can change size as a consequence of diet.Estrogen is produced by fat but this production is inhibited by AIs.Estrogen is also naturally washed out.In both diagrams, the lines ending with an arrow represent positive feedback whereas the lines ending with a bar denotes negative feedback.

Fig. 2 :
Fig. 2: Simulation results of the Eq.(2.1) for CD (left) and HFD (right) with the data points.Left axis corresponds to tumor size T (t) and fat volume F (t), right axis denotes estrogen concentration E(t) over time t.

Fig. 3 :Fig. 4 :
Fig. 3: Scenario Ia: Sum of the sensitive S and resistant R tumor subpopulations over time t for the constant treatment with different values of p associated with CD (left) and HFD (right).We mark the time point at which treatment is started with a dashed vertical line.

Fig. 5 :
Fig.5: Scenario Ia: Sum of the sensitive S and resistant R tumor subpopulations over time t for the optimal treatment with different values of p associated with CD (left) and HFD (right).The solid curve refers to the tumor size without treatment.We mark the time point at which treatment is started with a dashed vertical line.

Fig. 7 :
Fig. 7: Scenario Ia: Dynamics of model variables S, R, E and F over time t associated with CD (1st row ) and HFD (2nd row ).The solid curve refers to the tumor size without treatment, dotted curve corresponds to results for optimal treatment.

Fig. 8 :Fig. 9 :
Fig. 8: Scenario Ib: Sum of the sensitive S and resistant R tumor subpopulations over time t for the constant treatment with different values of p associated with CD (left) and HFD (right).We mark the time point at which treatment is started with a dashed vertical line.

Fig. 10 :Fig. 11 :
Fig.10: Scenario Ib: Sum of the sensitive S and resistant R tumor subpopulations over time t for the optimal treatment with different values of p associated with CD (left) and HFD (right).The solid curve refers to the tumor size without treatment.We mark the time point at which treatment is started with a dashed vertical line.

Fig. 12 :
Fig. 12: Scenario II: Sum of the sensitive S and resistant R tumor subpopulations over time t for the constant treatment with different values of p associated with CD (left) and HFD (right).We mark the time point at which treatment is started with a dashed vertical line.

Fig. 13 :
Fig.13: Scenario II: Left axis refers to the sum of the sensitive S and resistant R tumor subpopulations over time t for alternating treatment associated with CD (left) and HFD (right); right axis refers to treatment schedule.The dotte and solid curves refer to the tumor size with and without treatment, respectively.We mark the time point at which treatment is started with a dashed vertical line.

Fig. 14 :
Fig.14: Scenario II: Sum of the sensitive S and resistant R tumor subpopulations over time t for the optimal treatment with different values of p associated with CD (left) and HFD (right).The solid curve refers to the tumor size without treatment.We mark the time point at which treatment is started with a dashed vertical line.

Fig. 15 :
Fig.15: Scenario II: Optimal control function u over time t with ω R = ω S = ω U = 1.Dashed and dash-dotted curves refer to the optimal treatment schedules for CD and HFD, respectively.Solid line denotes the maximum treatment.We mark the time point at which treatment is started with dashed and dashdotted vertical lines for CD and HFD, respectively.

Fig. 16 :
Fig. 16: Scenario III: Dynamics of model variables S, R, E and F over time t associated with CD (1st row ) and HFD (2nd row ).

Fig. 17 :
Fig.17: Scenario III: Left axis refers to the sum of the sensitive S and resistant R tumor subpopulations over time t for alternating treatment associated with CD (left) and HFD (right); right axis refers to treatment schedule.The solid line refers to the tumor size without treatment.We mark the time point at which treatment is started with a dashed vertical line.

Fig. 18 :
Fig.18: Scenario III: Sum of the sensitive S and resistant R tumor subpopulations over time t for the optimal treatment with different values of p associated with CD (left) and HFD (right).The solid curve refers to the tumor size without treatment.We mark the time point at which treatment is started with a dashed vertical line.

Fig. 19 :
Fig.19: Scenario III: Optimal control function u over time t with ω R = ω S = ω U = 1.Dashed and dash-dotted curves refer to the optimal treatment schedules for CD and HFD, respectively.Solid line denotes the maximum treatment.We mark the time point at which treatment is started with dashed and dashdotted vertical lines for CD and HFD, respectively.

Fig. 20 :
Fig. 20: Scenario III: Dynamics of model variables S, R, E and F over time t associated with CD (1st row ) and HFD (2nd row ).The solid curve refers to the tumor size without treatment, dotted curve corresponds to results for optimal treatment.

3 √d
V is the number of layers with height d in the cube, where d is the diameter of an adipocyte.In addition, we compute 3 √

Table 1 :
Values of the parameters in the Eq.(2.1).