Complexity
Volume 2019 (2019), Article ID 9079104, 18 pages
https://doi.org/10.1155/2019/9079104

Research Article

Dynamic Modeling of the Angiogenic Switch and Its Inhibition by Bevacizumab

Dávid Csercsik1,2 and Levente Kovács¹

¹Physiological Controls Research Center, University Research, Innovation and Service Center, Óbuda University, Budapest, Hungary
²Pázmány Péter Catholic University, Faculty of Information Technology and Bionics, P.O. Box 278, 1444 Budapest, Hungary

Correspondence should be addressed to Dávid Csercsik; csercsik@itk.ppke.hu

Received 10 September 2018; Revised 12 November 2018; Accepted 12 December 2018; Published 3 January 2019

Guest Editor: Alejandro F. Villaverde

Copyright © 2019 Dávid Csercsik and Levente Kovács. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We formulate a dynamic model of vascular tumor growth, in which the interdependence of vascular dynamics with tumor volume is considered. The model describes the angiogenic switch; thus the inhibition of the vascularization process by antiangiogenic drugs may be taken into account explicitly. We validate the model against volume measurement data originating from experiments on mice and analyze the model behavior assuming different inputs corresponding to different therapies. Furthermore, we show that a simple extension of the model is capable of considering cytotoxic and antiangiogenic drugs as inputs simultaneously in qualitatively different ways.

1. Introduction

Neovascularization means the formation of new blood vessels. Angiogenesis, an important form of neovascularization, is characterized by hypoxia-driven sprouting of new capillaries from postcapillary venules. This mechanism plays an important role in many physiological (e.g., wound healing [1]) and pathological (e.g., macular degeneration [2]) processes. In the development of tumors, angiogenesis plays an exceptionally important role [3, 4]. In the beginning, when the tumorous cells form a small plaque, the tumor cells are well supported with metabolites by diffusion from the environment. However, as the size of the tumor increases, cells in the inside become insufficiently supported. Tumor-induced angiogenesis is the process of blood vessel formation, in which new vasculature is formed in order to support these insufficiently supported tumor cells.

Lately, much has been revealed about the details of tumor-induced angiogenesis and the underlying biochemical and biomechanical regulatory processes. These studies served as basis for the development of targeted molecular therapies [5]. The aim of these therapies is to inhibit tumor-related angiogenesis, thus cutting the tumor from metabolic support.

Bevacizumab (Avastin) is a pharmacotherapeutic antiangiogenic agent developed to withhold pathological angiogenesis [6] via the inhibition of the tumor angiogenic factor VEGF (vascular endothelial growth factor) [7]. VEGF may be considered as a representative member of the family of biochemical agents promoting angiogenesis, called tumor angiogenic factors (TAFs).

In [8], two different dosage protocols of bevacizumab were compared. In the case of the first protocol, experimental animals (mice) received one high dose according to the generally accepted medical principle, while in the case of the second protocol (quasi-continuous therapy), a much lower dose was delivered every day of the therapy. Results have shown that the quasi-continuous protocol was more effective, while the total injected amount of the drug was significantly less. As antiangiogenic agents are expensive, the total used drug amount is an important aspect to consider in therapeutic design. In addition, reduction of therapeutic doses is also desirable to minimize drug side effects. The result described in [8] underlines the importance of therapy optimization in the case of the application of antiangiogenic drugs and shows that computational methods may help to exploit the potentials of this approach.

On the other hand, based on the new data and paradigms brought to light in biological studies on angiogenesis, computational modeling of tumor-related vasculature development has became popular in the last decades, producing numerous computational models describing tumor growth and tumor-induced angiogenesis under different physiological circumstances (for a review on mathematical modeling of angiogenesis, see [9–12], while for a biological viewpoint on the mechanisms of physiological and tumor-related/pathological angiogenesis, one may refer to [13]). An important aspect of these modeling studies is to predict the effect of possible therapeutic approaches in cancer treatment [14–16].

A large part of the aforementioned models exhibit a quite high level of complexity, which implies that while they may be potentially appropriate for the comparison of different therapeutic approaches under the assumption of a given parameter set, it can be challenging to fit them to single patients. Furthermore, exact therapy optimization may be computationally infeasible if one relies on complex and spatially detailed models as [17, 18].

Feedback control [19, 20] may be an alternative to offline therapy optimization approaches. One benefit of closed-loop methods is that, assuming appropriate physiological signals, they may provide performance guarantees also in the presence of parametric uncertainties [21]. Other potential benefits of closed-loop treatments over protocol-based cancer therapies are discussed in [22]. In the case of diabetes, a similar biological control problem, such approaches have been successfully applied [23–27]. Control theoretic methods require, however, concentrated parameter models and ordinary differential (or difference) equation models with moderate complexity to perform well.

Recently, a simple dynamical model of tumor growth and the effect of the antiangiogenic drug bevacizumab has been published [28]. This model contains very few parameters and state variables and thus it is ideal for parameter estimation and controller design purposes. This model is based on three state variables, namely, the proliferating tumor volume, the necrotic tumor volume, and the concentration of the angiogenic inhibitor. Although this model provides a good fit for certain experimental data (see the data later in Section 2.2), it holds some flaws.

First, as it does not include the description of vasculature dynamics, its drawback is that it is unable to interpret advanced measurement data corresponding to tumor and vasculature evolution dynamics, potentially available in the foreseeable future. Recently, several imaging techniques have been described, which allow the reconstruction of vascular microstructures: Doppler optical frequency domain imaging [29] and functional photoacoustic microscopy [30] are used today already in in vivo setups to map vascular networks, while diffusible iodine-based contrast-enhanced computed tomography [31] may be used in terminal experimental animals. These methods could provide valuable data about vasculature dynamics in the near future, which may be used for the identification of the details of the angiogenic processes.

Second, minimal models not including vasculature dynamics as [28] are lacking the potential to describe the phenomenon of the angiogenic switch [32]. This hypothesis, formulated by Folkman, assumes that angiogenesis begins only at a certain stage of tumor development, more precisely at the time when the limited diffusion distance (which is about 0.1 mm [33]) makes the support of tumor cells inside the tumor with oxygen and metabolites no longer possible. According to the prediction of the minimal model [28], antiangiogenic drugs significantly affect the tumor growth also in the initial period. This contradicts with the consideration based on Folkman’s hypothesis that implies that in the initial period no insufficiently supported tumor cells are present, so no TAF synthesis is present; thus its effect can not be inhibited. In addition to the fact that minimal models like [28, 34] do not explicitly consider angiogenesis, the model [28] assumes a very simple constant rate of drug-independent necrosis, in which the proliferating cells turn into necrotic cells. In contrast, the process of necrosis strongly depends on the metabolic support and thus on the vascularization state of the tumor.

Third, paper [28] is based on the comparison of simulation results to measurement data originating from two scenarios. In these two scenarios, the antiangiogenic drug Bevacizumab was administered to experimental animals according to different protocols. In the first protocol, one 200 g bevacizumab dose was used for an 18-day therapy, while in the second (quasi-continuous) protocol one-tenth of the 200 g (20 g) dose was spread over 18 days; that is, 1.11 g bevacizumab was administered every day to the animals. Paper [8] is also based on results corresponding to these 2 administration protocols but also discusses results corresponding to therapy/drug-free case; namely, it states that mice that were treated with protocol 1 (one 200 g bevacizumab dose) did not have significantly smaller tumor volume than mice that did not receive therapy at all. In contrast to this result, if we compare the final tumor volumes resulting from the simulation of the model [28] in the case of protocol 1 and in the no-therapy case, we find that, in the case of protocol 1, the model prediction for the final tumor volume is 4741  (which is in good agreement with the experimental results), while in the no-therapy case, the model prediction for the final tumor volume is 37628 . This means that, according to the prediction of this model, the no-therapy final tumor size is almost 8 times larger than the protocol 1 case, making the validity of the model questionable in the no-therapy case, based on the results of [8]. We would underline that this does not question the validity of the model [28] in the case of Bevacizumab protocols and feedbacks similar to the ones discussed therein.

Another paper that aimed to formulate a control-oriented dynamical model was also recently published [35]. This model, which uses a bicompartmental approach, does describe the dynamics of vasculature but results in a 7-dimensional state space and 23 parameters, which is quite challenging size for control design. In addition, the model described in [35] was fit only to data corresponding to the first protocol in [8] (not both, like [28]).

According to the above preliminary results and considerations, our aim in this paper is as follows. We formulate a dynamic model, which includes the dynamics of the vasculature volume and describes the interplay between the tumor and vasculature volumes. To achieve this, we also include the dynamics of TAF in the model, which is produced by unsupported tumor and initiates the formation of new blood vessels from existing ones. This way our model will be capable of the interpretation of measurement corresponding to tumor and vasculature evolution dynamics. Furthermore, we aim to formulate a model, the predictions of which are also more acceptable in the therapy-free case.

2. Model Synthesis

In the following, we introduce the state-space variables of the model and interpret the state equations with the discussion of the model assumptions. Afterwards, the model parameters, their estimation, and the contextualization of the resulting parameters are described.

2.1. State Variables and Equations

The state-space variables of the model are summarized in Table 1 and their dynamics are described by the following equations:

Table 1: State variables of the model.

Notation
Variable
Dimension

V
Proliferating tumor volume

N
Necrotic tumor volume

B
Vasculature volume in the tumor

T
Concentration of TAF in the tumor
mg/ml

I
Inhibitor serum level
mg/ml

is the total volume of the tumor, the sum of the proliferating (living) volume and the necrotic volume: ; thus . Furthermore, denotes the unsupported living tumor volume ; thus . denotes the input, the injection of the antiangiogenic drug.

We assume a simple spherical tumor. This simplifying assumption (one-dimensional growth in other words) is widely used in the tumor modeling literature (see, e.g., [36–42]).

The auxiliary variable in the above equation is a key element of the model: it describes the actual ratio of the well-supported tumor cells in the tumor. is the function of the actual total volume and the vascularization ratio of the tumor (). can be computed as

where denotes the ideal density of vasculature in the tumor (notations with square brackets always refer to densities). This ideal density corresponds to the vasculature density, when all tumor cells are sufficiently supported; in other words, . In accordance with biological results [43–45], we assume that the vasculature is present in the living part of the tumor; thus the vasculature density is interpreted as density of blood vessels in the living (nonnecrotic) part of the tumor and is calculated as . Regarding the validity range of the model, we assume that the inequality holds at all times (in other words, we assume that the tumor is never fully vascularized).

The supported ratio of the tumor () is computed as

which may be viewed as a linear homotopy in between the constant 1 function and the function , which is also between zero and one. Thus also holds. The interpretation is that denotes a function that describes the ratio of tumor cells on the periphery of the tumor, which receive nutrients from outside of the tumor, so they are well supported. Naturally, this ratio depends on the actual total tumor volume . More precisely, in our terms, the periphery of the tumor is the outer shell of the sphere of our model, composed by tumor cells, which are closer to the surface of the tumor than the diffusion distance. According to [33], we assume the diffusion distance to be 150 . The ratio of periphery cells depends upon the radius of the tumor, which can be expressed from the tumor volume , according to the following simple derivation. Based on the tumor volume, we may derive the tumor radius (in mm) as

If , there is no tumor core. Assuming , the volume of the tumor core is

And thus the volume of the periphery is

and finally the function is derived as

If the tumor radius is below the diffusion distance, all cells are considered as periphery cells. The function describing the ratio of periphery cells in the range of 20-15000  is depicted in Figure 1 (according to [8], in the time of the first measurement, the volume of the tumors is about 50-60 ).

Figure 1: The function describing the ratio of periphery cells in the range of 20-15000  tumor volume.

Now let us return our focus to the auxiliary variable . The consideration behind the form of function (7) is that we assume two possible ways by which a tumor cell may get nutrient support. On one hand, if it is at the periphery of the tumor, it gets nutrients from the environment of the tumor via diffusion. On the other hand, if it is inside the tumor near a blood vessel, it also receives nutrient supply. The properties of function (7) reflect these considerations. We may see that if the tumor is composed only (or because of the approximation mostly) periphery cells or is almost fully vascularized (), the value of is approximately 1. In addition, at a fixed value of , it increases as the value of is increased, and at a fixed vale of , the value of is increased as increases (to put it simple, it is monotonically increasing in both variables).

Now, as we have discussed the interpretation of the auxiliary variables and the corresponding assumptions, we may return to the state equations. Equation (1) describes the tumor growth. The formula originates from the assumption that the well-supported part of the tumor volume () proliferates at the rate , while the unsupported volume necrotizes at the rate .

Equation (3) formalizes the dynamics of the vasculature, which may increase by two ways. On one hand, the term describes the internalization of new vasculature from the environment as the tumor grows, and on the other hand the term describes the formation of new blood vessels from existing ones in response to the TAF.

Equation (2) corresponds to necrotic volume. The importance of formally describing the process of necrosis lies in the fact that necrotized cells neither proliferate nor contribute to TAF production, so they must be distinguished from the general tumor volume.

Equations (4) and (5) describe the dynamics of TAF and the inhibitor concentration, respectively. The terms

of (4) describe that the production of TAF is proportional to the ratio of unsupported cells and takes place at rate , while it is cleared at the rate from the tumor. In the case of the inhibitor, the source is the injection (the input) , while based on [28], we assume that its clearance takes place according to Michaelis-Menten kinetics. The terms

and

in (4) and (5), respectively, correspond to the reaction in which the inhibitor binds the TAF molecule. Also, based on [28], the dynamics of the inhibition are considered assuming Michaelis-Menten kinetics with Michaelis-Menten constant (effective median dose).

The variable corresponds to the consideration that the concentrations of TAF () and the inhibitor () are not interpreted in the same volume (the compartment in the former case is the tumor; the volume in the latter case is the plasma). is the ratio of these volumes and can be computed as

where the value 1460 stands for the average blood volume of a mice in , in which the concentration of the inhibitor is interpreted [46].

2.2. Model Parameters

Some parameters of the model were taken from the literature (see Table 2), while the remaining parameters were estimated using experimental data originating from mice. Sápi et al. [8] carried out experiments, where C57Bl/6 mice with C38 colon adenocarcinoma were treated with bevacizumab using two different therapies.

(i)Therapy 1 (protocol-based treatment): five mice (mice C1-C5) were injected with 0.171 mg/ml bevacizumab at day 3 of the treatment (day 0 is considered as the day of the tumor implantation) (see Figure 2(a)).(ii)Therapy 2 (daily, quasi-continuous small amount administration): nine mice (mice E1-E9) received mg/ml injection of bevacizumab each day for 18 days from day 3 (see Figure 2(b)).

Table 2: Nominal model parameters () and their dimension, source, and 95% confidence interval (). PE stands for parameter estimation.

Notation
Dimension
Value
Source

1.4353
PE

-

PE

-

PE

PE

PE

PE

PE

PE

[48]
-

[49]
-

[28]
-

(a)

(b)

Figure 2: (a) The measured tumor volumes for mice C1-C5 that received therapy 1 and their average. (b) The measured tumor volumes for mice E1-E9 that received therapy 2 and their average.

The nominal parameter set of the model was determined using the average of measurements as reference and minimizing the mean square error of the deviance between the simulated and measured total volumes using the combination of particle swarm global optimization method [47] and the “nlinfit” function of MATLAB. During the simulations, the initial volume of the tumor was assumed to be equal to 5 , while the initial values of all other state variables were assumed to be 0. To capture the qualitatively different response of tumor growth to the two different protocols, both average curves (corresponding to protocols 1 and 2) were used simultaneously for the purpose of parameter estimation. While initial guess of the parameters was determined by the particle swarm global search method, the final values of the parameters were obtained using the “nlinfit” function. From the results of this function, the 95% confidence intervals () were determined using the function “nlparci.”

In order to potentially achieve a global optimum with the resulting parameter set, the estimation procedure was started from several initial coordinates in the parameter space. During parameter estimation, the average measurement results of both protocols were used to capture the qualitatively different response of the system to different inputs. Table 2 summarizes the model parameters. These parameters are to be referenced as nominal parameters in the further discussion and are denoted as .

In addition, to quantify parameter variance in the context of single trajectories, the model was fitted also for single growth curves as described in the Appendix.

Most estimated parameters of the model are hard to measure individually (in fact some of them are only interpreted inside the framework of this model) and no data are available on them which could serve as basis for comparison.

3. Results and Discussion

3.1. Fitting the Model to Measurement Data

In this subsection, we compare the model behavior and parameter values to measurement data. Figure 3 shows the fit of the model simulation output (total volume) to the experimental data sets that were used for the parameter estimation (therapy 1 and therapy 2). As it can be seen in the figures, the model sufficiently reproduces the measured growth trajectories, and the better fit is achieved in the case of therapy 2. For the sake of clarity, we note that since the average curve of measurement results was used for parameter estimation, the number of mice used in the experiments did not influence the objective function (in other words, this is not the reason for better fit in the case of therapy 2). To quantify the fit, we introduce the normalized squared deviation (NSD) as

where is the simulated total volume () and is the measured total volume. stands for the set of sample times, corresponding to days 3, 5, 7, 9, 11, 13, 15, 17, 19, and 21 in this case, while denotes the number of sample days. Normalization by the number of days is required, because we will use this measure later as well in the case of experimental data, where the number of sample points is lower. According to this measure, the representative values are and in the case of therapy 1 and therapy 2, respectively.

(a)

(b)

Figure 3: Measured and simulated tumor volumes in the case of therapy 1 (a) and therapy 2 (b).

The possible reason for the better fit in the case of the quasi-continuous therapy may be that the average of the measurements provides a more smooth, exponential-like curve in this case, which was possibly more easy to be achieved by the model.

The final volume of the tumor according to simulation results is 6670  in the case of therapy 1 and 3993  in the case of therapy 2. Considering the data described in Section 2.2, the final average measured volumes (corresponding to day 21) were 6604 and 3257  in the cases of therapy 1 and therapy 2, respectively.

3.2. Model Validation

3.2.1. Qualitative Validation

In this subsection, we analyze and compare the dynamic behavior of key model variables in the cases of no therapy, therapy 1, and therapy 2 defined in Section 2.2.

Figure 4 depicts the trajectory of and in the various cases.

(a)

(b)

Figure 4: (a) Ratio of the well-supported tumor cells () in the case of various therapies; (b) vascularization ratio () in the case of various therapies.

In general, it can be said that, in the beginning, when the tumor is small, due to the high ratio of periphery cells, (i.e., the ratio of the well-supported tumor cells) is high. However, as the tumor grows, the ratio of peripherial cells and so also decrease. The trajectories of the various cases differ only after day 6. The reason for this is that the process of angiogenesis becomes significant only if a large part of unsupported cells and TAF is already present, and the inhibition of the process becomes important only in this period. This is also shown in Figure 4(b), where it can be seen that the various cases regarding the vascularization ratio () differ only after day 6. Furthermore, it can be seen in Figure 4 that both and remain in the range as assumed during the model formulation. Figure 5 depicts the inhibitor concentration in the cases of therapy 1 and therapy 2 (in the case of no therapy, the inhibitor concentration is constant zero during the simulation). It can be seen in the figure that, in the case of therapy 1, the concentration is several orders of magnitude higher compared to therapy 2.

Figure 5: Concentration of the inhibitor () in the case of various therapies on two different scales.

Figure 6 depicts the concentration of TAF () and the term of (3) (corresponding ot TAF-dependent vascularization) in response to the various levels of the inhibitor in the case of different therapies.

(a)

(b)

Figure 6: (a) Concentration of TAF () in the case of various therapies; (b) TAF-dependent vascularization (TdV) in the case of different therapies.

In this figure (Figure 6(a)), it can be seen that the concentration of TAF () around days 3-8 is almost the same in the case of therapy 1 and therapy 2, although the concentration of the inhibitor () differs on the orders of magnitude. The explanation for this is the effective median dose () of the inhibitor, which defines the inhibitor level at which the effect of the drug saturates. We have to note that this parameter was not subject to estimation but it was taken from the article [48]. Figure 6(b) depicts that from the time when the large bolus has been cleared in the case of therapy 1 (around day 13), the TAF-dependent vascularization (TdV) becomes significantly different in the case of the two therapies.

Since the volume trajectories of the model follow an exponential-like fashion in all cases, they are maybe not so informative as, for example, the plot of , but for the sake of completeness, they are depicted in Figure 7.

(a)

(b)

(c)

(d)

Figure 7: Volume trajectories of the model (living volume , unsupported volume , necrotic volume , and total volume , corresponding to subfigures (a), (b), (c), and (d), resp.) in the case of various therapies.

3.2.2. Validation against Measurement Data with No Therapy

In phase I of the experiments described in [8], the experimental animals received no treatment. In this case, the length of the experiment was only 19 days (compared to 21 days in all other experiments discussed before and depicted in Figure 2). The measurement data originating from this experiment was used to validate the model. Figure 8 shows the results of the model validation.

Figure 8: Measured and simulated tumor volumes in the case of no therapy.

The measure of the fit introduced in (16) gives a value of in this case (the lower index corresponds to “no therapy”). If one compares this value to the NSD values corresponding to the experimental data to which the model was fitted (see Section 3.1), or examines the figures, it can be seen that the error of the fit is one order of magnitude higher in this case (since the no-therapy case was not considered during parameter estimation). On the other hand, taking into account the significant variation among experimental animals as well, the the model output in the case of no therapy provides an acceptable fit with measurement data.

3.2.3. Validation regarding Final Tumor Volume Values in Various Cases

As in control applications, for which the current model is primary proposed, the usual aim is to minimize the final volume of the tumor (under, e.g., constraints corresponding to the total applied drug quantity); it is important to compare the final tumor sizes. As the final day of the experiments was different in the no-therapy case and therapies 1 and 2, Tables 3 and 4 summarize the simulated and measured tumor volumes on days 19 and 21.

Table 3: Simulated and measured average total volumes on day 21 ().

Therapy

No therapy
1.03
No data

Therapy 1
6.67
6.60

Therapy 2
3.99
3.26

Table 4: Simulated and measured average total volumes on day 19 ().

Therapy

No therapy
4.68
6.15

Therapy 1
3.10
4.44

Therapy 2
2.13
2.03

For the sake of comparison to previous literature results, let us note that the simulated volume on day 19 assuming no therapy is in the case of the model described in [28] (this value has been obtained by the reproduction of the model described in [28]). Comparing the differences of the values, we may say that the validity of the proposed model compared to [28] is significantly better regarding the final volume in the no-therapy case.

3.3. Model Identifiability

In the following two subsections, we present some results related to the parameter sensitivity and structural identifiability of the proposed model.

3.3.1. Parameter Sensitivity of the Model

In this subsection, we analyze the parameter sensitivity of the model for the estimated parameters. The sensitivity analysis is an important tool to characterize how the model parameters affect the simulation output. The presence of very large differences in the sensitivities of parameters may point to identifiability problems.

In order to formalize this analysis, we define the sensitivity measure detailed in the following equation:

In (17), denotes the nominal parameter vector detailed in Table 2, while denotes a perturbed parameter vector. and stand for the nominal output () of the model and for the output in the case of the perturbed parameter vector, respectively.

No-Therapy Case. Table 5 summarizes the results in the no-therapy case. Each row of this table corresponds to a perturbed parameter vector, in which only one element differs from the nominal vector (by 20, 10, or 5%).

Table 5: Sensitivities () of the model to the changes of parameters , , , , , , , , , and in the no-therapy case. One unit is .

-20%
-10%
-5%
+5%
+10%
+20%

0.1321
0.0240
0.0052
0.0040
0.0143
0.0462

0.4715
0.1601
0.0467
0.0640
0.2996
1.6456

0.0653
0.0173
0.0044
0.0047
0.0193
0.0821

0.0072
0.0017
0.0004
0.0004
0.0016
0.0061

0.2287
0.0769
0.0224
0.0309
0.1454
0.8147

0.2287
0.0769
0.0224
0.0309
0.1454
0.8148

0.5561
0.0868
0.0176
0.0120
0.0405
0.1185

0
0
0
0
0
0

First, it is conspicuous that the sensitivity to the parameter is zero in this case. The explanation for this is that this parameter corresponds to the effect of the angiogenic inhibitor (), the concentration of which is constantly 0 in this case. In other words, there is no drug effect in this case, the dynamics of which are affected by this parameter.

Second, we can see in Table 5 that the sensitivity for is equal to the sensitivity of in the no-therapy case. The reason for this is that while corresponds to the synthesis rate of TAF, corresponds to the rate of TAF-dependent vascularization. If no inhibitor is present, there is no difference between an increase in the TAF concentration and a more efficient TAF-driven vascularization (see (3) and (4)). Let us note that if the inhibitor is present, not the whole portion of the synthetized TAF takes part in the blood vessel formation process (since some molecules are binding to the inhibitor); thus the situation is different if input (thus nonzero drug concentration) is also present.

Third, the model is the most sensitive to the parameter . This is not surprising, as directly affects the dynamics of (see (1)) as a proportional term, so its effect in , which is directly present in the output, is exponential. Apart form , , and , the model shows the least sensitivity to , the rate of necrosis.

Therapy 1. Table 6 shows that the model sensitivity for , the parameter corresponding to the effect of the inhibitor, is relatively low.

Table 6: Sensitivities () of the model to the changes of parameters , , , , , , , , , and in the case of therapy 1. One unit is .

-20%
-10%
-5%
+5%
+10%
+20%

0.2926
0.0525
0.0114
0.0086
0.0305
0.0975

0.9383
0.3259
0.0962
0.1353
0.6422
3.6265

0.1375
0.0368
0.0095
0.0102
0.0420
0.1807

0.0181
0.0043
0.0011
0.0010
0.0039
0.0149

0.4498
0.1518
0.0442
0.0611
0.2889
1.6242

0.4594
0.1559
0.0456
0.0634
0.3011
1.7086

0.9364
0.1493
0.0305
0.0213
0.0720
0.2128

0.0100
0.0023
0.0005
0.0005
0.0019
0.0069

Table 7 shows that, in the case of therapy 2, an extreme high sensitivity is experienced in the case of the increase of parameter . Moreover, the sensitivity to is also significantly decreased. As the injections and thus the concentrations of the inhibitor () are by 2 orders of magnitude lower compared to therapy one (see Figure 5), it is plausible that the exact value of its clearance parameter () has a significantly less effect on the dynamics of the model compared to therapy 1, where a large dose is applied.

Table 7: Sensitivities () of the model to the changes of parameters , , , , , , , , , and in the case of therapy 2. One unit is .

-20%
-10%
-5%
+5%
+10%
+20%

0.1304
0.0230
0.0050
0.0038
0.0131
0.0415

0.4012
0.1422
0.0427
0.0620
0.3015
1.7891

0.0582
0.0157
0.0041
0.0044
0.0185
0.0798

0.0090
0.0021
0.0005
0.0005
0.0019
0.0072

0.1232
0.0402
0.0116
0.0155
0.0722
0.3971

0.1332
0.0448
0.0131
0.0179
0.0859
0.4902

0.1690
0.0288
0.0061
0.0045
0.0154
0.0473

0.0135
0.0030
0.0007
0.0006
0.0023
0.0081

Apart from this, the results are similar to the case of therapy 1.

Altogether, based on the results of the sensitivity analysis, it can be said that further experiments focusing solely on pharmacokinetics of the applied drugs are desirable to estimate the parameter decoupled from other model parameters.

3.3.2. Structural Identifiability

Structural identifiability properties of a system describe whether there is a theoretical possibility for the unique determination of system parameters from appropriate input-output measurements or not. It is important to emphasize that identifiability is a property of the model structure. Basic early references for studying identifiability of dynamical systems are [50, 51]. Since the introduction of this concept, multiple approaches have been proposed for the analysis of structural identifiability of various nonlinear system classes, for example, polynomial systems [52] or rational function state-space models [53]. A critical comparison of methods for identifiability analysis is described in [54].

First, let us note that as our model uses the variable describing vasculature density, our model falls into the class of rational function state-space models.

Second, let us consider the factors that make the structural identifiability analysis challenging in our case. Structural identifiability methods usually rely on iterative computation of (Lie-) derivatives of the output (see, e.g., [55], on which the software used later is based).

For identifiability analysis, let us consider a reduced version of the proposed model, which assumes no input (no antiangiogenic drug is present). The simplified form of the model is described by (18)-(21). As we will see, this submodel already poses a challenge regarding identifiability due to the complexity of the resulting equations.

In this case, if only the total volume may be measured (as in the case of our measurements used for the parameter estimation),

If we consider further derivatives,

where

It is easy to se that, in , on one hand, the derivative of the function appears, and, on the other hand, the derivative of is also present. These are long and complicated expressions, the higher-order derivatives of which are needed in the further steps.

Based on the above considerations, for the structural identifiability analysis, we use the freely available GenSSI [56, 57] software, which is able to handle complex expressions with the help of computer-algebra methods. GenSSI implements iteratively the generating series method, as presented in [58], with the help of identifiability tableaus, as described in [55].

According to the results of this software, the parameter of the model is structurally globally identifiable, but neither positive (structural global/local identifiability) nor negative (structural nonidentifiability) results are obtained for other parameters. This result is based on 7th-order Lie-derivatives, which has been proven to be the computational limit in our case.

Nevertheless, let us discuss this topic a bit further from the point of view of possible future measurements with regard to the proposed model. In the recent years, multiple imaging techniques have been developed, which allow the 3D reconstruction of vascular microstructures: Doppler optical frequency domain imaging [29] and functional photoacoustic microscopy [30] are used today already in in vivo setups to map vascular networks, while diffusible iodine-based contrast-enhanced computed tomography [31] may be used in terminal experimental animals.

If these methods will be applicable in the case of animals used in the experiments, total tumor volume will be accessible during measurements together with the total vasculature volume . Interpreted for our case, this will mean that we will have two observables and . If we rerun the structural identifiability analysis with this new model output, we get the result that , , , and are structurally globally identifiable (and no result is obtained for other parameters, similar to the previous case, so they might be or might be not identifiable). In this case, the maximum order of the Lie-derivatives, for which the computation was feasible, was the 6th order.

The complete identifiability tableaus of the reduced model are depicted in Figure 9. For the interpretation of these tableaus, see [55] or the GenSSI UserGuide [56].

(a)

(b)

Figure 9: Complete identifiability tableaus of the reduced model in the case when the output is (a) and in the case when the output is (b).

Based on the above, it may be suspected that the model will have beneficial properties shall it be fitted for measurements planned to be carried out in the foreseeable future.

3.4. Extension of the Model in order to Account for Combined Therapy

In the clinical practice, antiangiogenic drugs are often used together with conventional cytotoxic substances. In this setup, while the cytotoxic agent enhances the degeneration/necrosis of tumor cells, the antiangiogenic drugs are resposible for cutting the tumor from metabolic support via the inhibition of angiogenesis. Several results have been published recently corresponding to these combined therapies [59, 60].

Models with predictive power regarding the efficiency of combined therapies and model-based optimization of such treatments are not prevalent in literature. Some initial results on the optimization of combined therapies are described in [61], using the model of [35].

As the proposed model is taking into account vasculature and tumor cell dynamics in a differentiated way, it is able to distinguish between qualitatively different inputs related to different therapeutic agents. As a consequence, the proposed model may be easily extended to consider not only angiogenic drugs but also cytotoxic drugs. Let us consider the following modified state-space model described in the following equations:

First, the new equation (31) describes the time evolution of the cytotoxic drug, the injection of which is described by the term . To clarify notations, the injection of the angiogenic inhibitor is denoted by in this case. The term describes the clearance of the cytotoxic drug; the parameter denotes its clearance rate. In this case, we assume a simple clearance (no saturation dynamics). The reason for this is on one hand that this approach requires less parameters, and on the other hand as long as the exact identity of the cytotoxic drug is unknown, the dynamical features of its clearance can not be precisely determined (of course the clearance dynamics may be later refined).

The effect of the cytotoxic drug is modeled in this case as an enzymatic reaction, in which the cytotoxic drug acts as an enzyme, turning living cells to necrotic cells. This mechanism is described by the term

in (26) and by the complementary term

in (27). and are new parameters describing the efficiency of the cytotoxic drug in enzymatic context assuming Michaelis-Menten kinetics.

This way the effects of the two drugs are considered in qualitatively different ways in the model. While the antiangiogenic drug acts explicitly on the formation of new blood vessels by binding to TAF and thus inhibiting angiogenesis, the cytotoxic drug acts as an enzyme, driving living tumor cells to necrosis, independent of the actual vascular state of the tumor.

4. Conclusions and Future Work

In this article, we formulated a dynamic model of vascular tumor growth, which accounts for the vasculature and TAF concentration development of the tumor and thus is able to reproduce the phenomenon of the angiogenic switch. We validated the model against volume measurement data originating from experiments on mice and found that the model provides a good fit for tumor volume data in both cases of the two analyzed therapies. The extension of the model described in Section 3.4 makes the model capable of accounting for qualitatively different effect of antiangiogenic and cytotoxic drugs.

When comparing the proposed model to literature results, we may state the following. Regarding the model in [28] (considering the extended model described therein), the proposed approach uses more state variables (5 vs 3) and holds more parameters (12 versus 8) but describes vasculature dynamics as well. This feature will allow us to fit the model to dynamical vasculature data, hopefully available in the foreseeable future, and thus get a more precise dynamical representation of angiogenesis-dependent tumor growth and its inhibition. Furthermore, the validity of the proposed model compared to [28] seems better regarding the no-therapy case.

Comparing the model described in the current article to [35], we see that although the model described in [35] accounts for vasculature dynamics as well, it uses more state variables (7 versus 5) and significantly more parameters (22 versus 12). Furthermore, the model in [35] was fitted only for measurement data originating from protocol 1, while the proposed model has been validated against both protocol 1 and protocol 2. Similar to the model described in [35], the model proposed in the current article also allows for the analysis of combined therapies, as done in the article in [61].

Regarding future work, in the framework of the project Tamed Cancer (ERC grant agreement number 679681), animal experiments (mice) aiming to characterize the vasculature development during tumor growth are planned in the near future. These experiments will provide reference data for both vasculature volumes and tumor volumes, so we will be able to fit the model in either dimension against experimental data. This will allow further validation, refinement, or recalibration of the model.

Experiments regarding the efficiency of various combined therapies are also expected in the future, which will serve as reference scenarios regarding the identification of the extended model described in Section 3.4.

Once the model is identified and validated from multiple aspects, studies on therapy optimization in open-loop and closed-loop setup will take place.

Appendix

In this appendix, we detail the fitting of the model to the individual trajectories corresponding to single mice in the case of either therapy 1 or therapy 2. was assumed in all cases. Figures 10(a)–10(e) depict the fit of the model output () to the individual measured growth trajectories in the case of therapy 1, while Figures 11(a)–11(i) depict the fit of the model output to the individual measured growth trajectories in the case of therapy 2.

(a)

(b)

(c)

(d)

(e)

Figure 10: Fitting the model to individual trajectories in the case of therapy 1: simulated output (with the parameters obtained from fitting the model to the specific trajectory (), simulated output with nominal parameters (), and measured output.)

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

Figure 11: Fitting the model to individual trajectories in the case of therapy 2: simulated output (with the parameters obtained from fitting the model to the specific trajectory (), simulated output with nominal parameters (), and measured output.)

In every case, the simulated output assuming the nominal parameters detailed in Table 2 is also depicted to serve as basis of comparison. Tables 8 and 9 summarize the parameter values resulting in the case of fitting to individual trajectories in the case of therapy 1 and therapy 2, respectively.

Table 8: Model parameter values resulting from fitting to individual trajectories in the case of therapy 1 (mice -).

Parameter

Scale

0.0336
0.0127
0.0289
0.0374
0.0246

1.4336
1.4454
1.55
1.435
1.4154

1.4403
1.9642
1.4604
2.2084
1.2167

0.0838
0.1684
0.0737
0.0739
0.0791

4.4337
5.5455
10.2252
4.4032
4.4909

1.8021
0.5612
0.7111
1.9291
1.4795

0.2762
0.3917
0.2517
0.1173
0.3162

0.1995
0.1637
0.1701
0.2114
0.1369

Table 9: Model parameter values resulting from fitting to individual trajectories in the case of therapy 2 (mice -).

Par.

Scale

0.0322
0.0393
0.0331
0.0184
0.039
0.0179
0.0422
0.0223
0.0331

1.622
0.9568
1.436
1.433
1.408
1.442
1.634
1.609
1.432

1.463
2.01
1.373
2.026
1.195
2.004
1.173
0.9805
0.6332

0.0933
0.0552
0.1056
0.0567
0.0618
0.1236
0.0933
0.0425
0.1056

17.34
18.67
18.08
18.66
18.66
18.08
18.08
18.65
18.08

0.444
0.6039
0.4612
0.4512
0.923
0.1138
0.6959
0.335
0.9693

0.2834
0.2657
0.4071
0.3201
0.3063
0.4105
0.3123
0.3311
0.4066

0.2449
0.1011
0.2605
0.2546
0.259
0.2496
0.1749
0.2643
0.2603

Table 10 holds the standard deviation (STD) values of the estimated parameters obtained by fitting the model to individual trajectories. Let us emphasize that, in the case of these estimates, only one of the protocols was considered for the fitting, namely, the one from which the actual trajectory originates. Naturally, it is a significantly harder task to find a parameter set for the model which describes the response to both protocols simultaneously (as in Table 2). Parameters obtained by fitting to single volume trajectories are potentially unable to appropriately describe the response to the two protocols at the same time. The relatively high STD values are not surprising in the light of the significant differences among individual tumor growth trajectories.

Table 10: Standard deviation (STD) of the estimated parameters, regarding fitting to individual trajectories compared to their nominal value in %.

Parameter
STD (%)

10.96

23.40

39.60

33.33

52.64

47.20

31.40

32.57

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant Agreement no. 679681).

Supplementary Materials

simulate_MMAGS.m is the file that performs the simulation of the model in MATLAB; it uses the following files: inj_fnc_prot_1_discrete.m, description of the injection function in the case of therapy 1; inj_fnc_prot_2_discrete.m, description of the injection function in the case of therapy 2; f_gamma.m, implementation of the function ∖gamma detailed in the article (see (7)). sim_tumor_dynamics_MMAGS_2_discrete.m is one central file that calls simulate_MMAGS.m; it has two important input parameters, flag_plot_details, may be 0 or 1, if it is one, more figures are generated, protocol - can be 0, 1 or 2 depending on the actual simulated protocol (therapy); it depicts the simulation results and compares them with experimental data (Figure 3 in the manuscript). Experimental data is defined inside the file. Corresponding to protocol 1, the variable “adatsor” holds the vector of the experimental data. Corresponding to protocol 2, the variable “adatsor_2” holds the vector of the experimental data. In both cases, the data corresponding to various mice are integrated in one vector (see its decomposition later in the file. In the case of protocol 1, the variable sote_data holds the decomposed version; in the case of protocol 2, the variables “Eger1” and so forth hold the decomposed data). sim_compare_discrete_1.m performs the simulations in the case of no therapy, protocol 1, and protocol 2 and depicts the results corresponding to the 3 different cases in figures (Figures 4–7 in the manuscript). (Supplementary Materials)

References

P. C. Lee, A. N. Salyapongse, G. A. Bragdon et al., “Impaired wound healing and angiogenesis in eNOS-deficient mice,” American Journal of Physiology-Heart and Circulatory Physiology, vol. 277, no. 4, pp. H1600–H1608, 1999.

E. W. M. Ng and A. P. Adamis, “Targeting angiogenesis, the underlying disorder in neovascular age-related macular degeneration,” Canadian Journal of Ophthalmology, vol. 40, no. 3, pp. 352–368, 2005.

J. Folkman, “The role of angiogenesis in tumor growth,” Seminars in Cancer Biology, vol. 3, no. 2, pp. 65–71, 1992.

J. Folkman, “Role of angiogenesis in tumor growth and metastasis,” Seminars in Oncology, vol. 29, 6, supplement 6, pp. 15–18, 2002.

T. A. Baudino, “Targeted cancer therapy: the next generation of cancer treatment,” Current Drug Discovery Technologies, vol. 12, no. 1, pp. 3–20, 2015.

N. Ferrara, K. J. Hillan, and W. Novotny, “Bevacizumab (Avastin), a humanized anti-VEGF monoclonal antibody for cancer therapy,” Biochemical and Biophysical Research Communications, vol. 333, no. 2, pp. 328–335, 2005.

N. Ferrara, H. Gerber, and J. LeCouter, “The biology of VEGF and its receptors,” Nature Medicine, vol. 9, no. 6, pp. 669–676, 2003.

J. Sápi, L. Kovács, D. A. Drexler, P. Kocsis, D. Gajári, and Z. Sápi, “Tumor volume estimation and quasi- continuous administration for most effective bevacizumab therapy,” PLoS ONE, vol. 10, no. 11, Article ID e0142190, pp. 1–20, 2015.

M. A. J. Chaplain, S. R. McDougall, and A. R. A. Anderson, “Mathematical modeling of tumor-induced angiogenesis,” Annual Review of Biomedical Engineering, vol. 8, pp. 233–257, 2006.

S. M. Peirce, “Computational and mathematical modeling of angiogenesis,” Microcirculation, vol. 15, no. 8, pp. 739–751, 2008.

M. Scianna, C. G. Bell, and L. Preziosi, “A review of mathematical models for the formation of vascular networks,” Journal of Theoretical Biology, vol. 333, pp. 174–209, 2013.

H. Rieger and M. Welter, “Integrative models of vascular remodeling during tumor growth,” Wiley Interdisciplinary Reviews: Systems Biology and Medicine, vol. 7, no. 3, pp. 113–129, 2015.

M. Papetti and I. M. Herman, “Mechanisms of normal and tumor-derived angiogenesis,” American Journal of Physiology-Cell Physiology, vol. 282, no. 5, pp. C947–C970, 2002.

H. A. Levine, B. D. Sleeman, and M. Nilsen-Hamilton, “A mathematical model for the roles of pericytes and macrophages in the initiation of angiogenesis. I. The role of protease inhibitors in preventing angiogenesis,” Mathematical Biosciences, vol. 168, no. 1, pp. 77–115, 2000.

L. Arakelyan, V. Vainstein, and Z. Agur, “A computer algorithm describing the process of vessel formation and maturation, and its use for predicting the effects of anti-angiogenic and anti-maturation therapy on vascular tumor growth,” Angiogenesis, vol. 5, no. 3, pp. 203–214, 2002.

J. Poleszczuk, P. Hahnfeldt, and H. Enderling, “Therapeutic implications from sensitivity analysis of tumor angiogenesis models,” PLoS ONE, vol. 10, no. 3, 2015.

A. L. Bauer, T. L. Jackson, and Y. Jiang, “A cell-based model exhibiting branching and anastomosis during tumor-induced angiogenesis,” Biophysical Journal, vol. 92, no. 9, pp. 3105–3121, 2007.

J. L. Gevertz and S. Torquato, “Modeling the effects of vasculature evolution on early brain tumor growth,” Journal of Theoretical Biology, vol. 243, no. 4, pp. 517–531, 2006.

G. F. Franklin, D. J. Powell, and A. Emami-Naeini, Feedback control of dynamic systems, vol. 3, Addison-Wesley, Reading, MA, USA, 1994.

J. C. Doyle, B. A. Francis, and A. R. Tannenbaum, Feedback control theory, Courier Corporation, 2013.

Z. Kemin and J. C. Doyle, Essentials of Robust Control, vol. 104, Prentice Hall, Upper Saddle River, NJ, USA, 1998.

J. Sápi, D. A. Drexler, and L. Kovács, “Potential benefits of discrete-time controller-based treatments over protocol-based cancer therapies,” Acta Polytechnica Hungarica, vol. 14, no. 1, pp. 11–23, 2017.

J. Bondia, S. Romero-Vivo, B. Ricarte, and J. L. Diez, “Insulin estimation and prediction: a review of the estimation and prediction of subcutaneous insulin pharmacokinetics in closed-loop glucose control,” IEEE Control Systems Magazine, vol. 38, no. 1, pp. 47–66, 2018.

F. J. Doyle, L. M. Huyett, J. B. Lee, H. C. Zisser, and E. Dassau, “Closed-loop artificial pancreas systems: Engineering the algorithms,” Diabetes Care, vol. 37, no. 5, pp. 1191–1197, 2014.

R. A. Harvey, E. Dassau, W. C. Bevier et al., “Clinical evaluation of an automated artificial pancreas using zone-model predictive control and health monitoring system,” Diabetes Technology & Therapeutics, vol. 16, no. 6, pp. 348–357, 2014.

L. Kovács, “Linear parameter varying (LPV) based robust control of type-I diabetes driven for real patient data,” Knowledge-Based Systems, vol. 122, pp. 199–213, 2017.

P. H. Colmegna, R. S. Sanchez-Pena, R. Gondhalekar, E. Dassau, and F. J. Doyle, “Switched LPV Glucose Control in Type 1 Diabetes,” IEEE Transactions on Biomedical Engineering, vol. 63, no. 6, pp. 1192–1200, 2016.

D. A. Drexler, J. Sápi, and L. Kovács, “Modeling of Tumor Growth Incorporating the Effects of Necrosis and the Effect of Bevacizumab,” Complexity, vol. 2017, Article ID 5985031, 10 pages, 2017.

B. J. Vakoc, R. M. Lanning, J. A. Tyrrell et al., “Three-dimensional microscopy of the tumor microenvironment in vivo using optical frequency domain imaging,” Nature Medicine, vol. 15, no. 10, pp. 1219–1223, 2009.

H. F. Zhang, K. Maslov, G. Stoica, and L. V. Wang, “Functional photoacoustic microscopy for high-resolution and noninvasive in vivo imaging,” Nature Biotechnology, vol. 24, no. 7, pp. 848–851, 2006.

P. M. Gignac, N. J. Kley, J. A. Clarke et al., “Diffusible iodine-based contrast-enhanced computed tomography (diceCT): An emerging tool for rapid, high-resolution, 3-D imaging of metazoan soft tissues,” Journal of Anatomy, vol. 228, no. 6, pp. 889–909, 2016.

D. Hanahan and J. Folkman, “Patterns and emerging mechanisms of the angiogenic switch during tumorigenesis,” Cell, vol. 86, no. 3, pp. 353–364, 1996.

S. Redline and N. A. Berger, Impact of Sleep and Sleep Disturbances on Obesity and Cancer, Springer, New York, NY, USA, 2014.

D. A. Drexler, J. Sápi, and L. Kovács, “A minimal model of tumor growth with angiogenic inhibition using bevacizumab,” in Proceedings of the 15th IEEE International Symposium on Applied Machine Intelligence and Informatics, SAMI 2017, pp. 185–190, Slovakia, January 2017.

D. Csercsik, J. Sapi, T. Gonczy, and L. Kovacs, “Bi-compartmental modelling of tumor and supporting vasculature growth dynamics for cancer treatment optimization purpose,” in Proceedings of the 2017 IEEE 56th Annual Conference on Decision and Control (CDC), pp. 4698–4702, Melbourne, VIC, December 2017.

J. Landry, J. P. Freyer, and R. M. Sutherland, “A model for the growth of multicellular spheroids,” Cell Proliferation, vol. 15, no. 6, pp. 585–594, 1982.

S. A. Maggelakis and J. A. Adam, “Mathematical model of prevascular growth of a spherical carcinoma,” Mathematical and Computer Modelling, vol. 13, no. 5, pp. 23–38, 1990.

G. W. Swan, “The diffusion of an inhibitor in a spherical tumor,” Mathematical Biosciences, vol. 108, no. 1, pp. 75–79, 1992.

M. A. J. Chaplain, D. L. Benson, and P. K. Maini, “Nonlinear diffusion of a growth inhibitory factor in multicell spheroids,” Mathematical Biosciences, vol. 121, no. 1, pp. 1–13, 1994.

M. A. J. Chaplain, “Avascular growth, angiogenesis and vascular growth in solid tumours: The mathematical modelling of the stages of tumour development,” Mathematical and Computer Modelling, vol. 23, no. 6, pp. 47–87, 1996.

S. A. Maggelakis, “The effects of tumor angiogenesis factor (TAF) and tumor inhibitor factors (TIFs) on tumor vascularization: A mathematical model,” Mathematical and Computer Modelling, vol. 23, no. 6, pp. 121–133, 1996.

H. M. Byrne, J. R. King, D. L. McElwain, and L. Preziosi, “A two-phase model of solid tumour growth,” Applied Mathematics Letters, vol. 16, no. 4, pp. 567–573, 2003.

R. H. Thomlinson and L. H. Gray, “The histological structure of some human lung cancers and the possible implications for radiotherapy,” British Journal of Cancer, vol. 9, no. 4, pp. 539–549, 1955.

K. I. E. M. Wÿffels, J. H. A. M. Kaanders, P. F. J. W. Rijken et al., “Vascular architecture and hypoxic profiles in human head and neck squamous cell carcinomas,” British Journal of Cancer, vol. 83, no. 5, pp. 674–683, 2000.

J. Forster, W. Harriss-Phillips, M. Douglass, and E. Bezak, “A review of the development of tumor vasculature and its effects on the tumor microenvironment,” Hypoxia, vol. Volume 5, pp. 21–32, 2017.

J. J. Heatley, “Cardiovascular Anatomy, Physiology, and Disease of Rodents and Small Exotic Mammals,” Veterinary Clinics of North America - Exotic Animal Practice, vol. 12, no. 1, pp. 99–113, 2009.

A. I. Vaz and L. N. Vicente, “A particle swarm pattern search method for bound constrained global optimization,” Journal of Global Optimization, vol. 39, no. 2, pp. 197–219, 2007.

H.-P. Gerber and N. Ferrara, “Pharmacology and pharmacodynamics of bevacizumab as monotherapy or in combination with cytotoxic therapy in preclinical studies,” Cancer Research, vol. 65, no. 3, pp. 671–680, 2005.

European Medicines Agency. Scientific discussion. 2005. http://www.ema.europa.eu/docs/en_GB/document_library/EPAR_-_Scientific_Discussion/human/000582/WC500029262.pdf.

E. Walter, Identifiability of state space models, Springer, Berlin, Germany, 1982.

E. Walter, Identifiability of Parametric Models, Elsevier, 2014.

S. Vajda, “Identifiability of polynomial systems: structural and numerical aspects,” Identifiability of Parametric Models, vol. 4, pp. 42–49, 1987.

G. Margaria, E. Riccomagno, M. . Chappell, and H. P. Wynn, “Differential algebra methods for the study of the structural identifiability of rational function state-space models in the biosciences,” Mathematical Biosciences, vol. 174, no. 1, pp. 1–26, 2001.

O.-T. Chis, J. R. Banga, and E. Balsa-Canto, “Structural identifiability of systems biology models: A critical comparison of methods,” PLoS ONE, vol. 6, no. 11, p. e27755, 2011.

E. Balsa-Canto, A. Antonio Alonso, and R. Julio Banga, “An iterative identification procedure for dynamic modeling of biochemical networks,” BMC systems Biology, vol. 4, no. 1, p. 11, 2010.

O. Chiş, J. R. Banga, and E. Balsa-Canto, “GenSSI: A software toolbox for structural identifiability analysis of biological models,” Bioinformatics, vol. 27, no. 18, Article ID btr431, pp. 2610-2611, 2011.

T. S. Ligon, F. Fröhlich, O. T. Chiş et al., “GenSSI 2.0: multi-experiment structural identifiability analysis of SBML models,” Bioinformatics, vol. 34, no. 8, pp. 1421–1423, 2018.

E. Walter and Y. Lecourtier, “Unidentifiable compartmental models: what to do?” Mathematical Biosciences, vol. 56, no. 1-2, pp. 1–25, 1981.

M. Scalerandi, B. C. Sansone, and C. A. Condat, “Diffusion with evolving sources and competing sinks: development of angiogenesis,” Physical review. E, Statistical, nonlinear, and soft matter physics, vol. 65, no. 1, Article ID 011902, p. 011902, 2002.

R. K. Jain, “Normalizing tumor vasculature with anti-angiogenic therapy: A new paradigm for combination therapy,” Nature Medicine, vol. 7, no. 9, pp. 987–989, 2001.

D. Csercsik, J. Sapi, and L. Kovacs, “Model-Based Simulation and Comparison of Open-Loop and Closed-Loop Combined Therapies for Tumor Treatment,” in Proceedings of the 2018 IEEE Conference on Control Technology and Applications (CCTA), pp. 1383–1388, Copenhagen, August 2018.