CANCER DISEASE
INTEGRATIVE MODELLING APPROACHES
Helen M. Byrne, Markus R. Owen∗,
Centre for Mathematical Medicine,
School of Mathematical Sciences,
University of Nottingham,
Nottingham NG7 2RD,
England,
Tomas Alarcon† and Philip K. Maini,
Centre for Mathematical Biology,
Mathematical Institute,
24-29 St Giles’,
Oxford OX1 3LB,
England.
Centre for Mathematical Medicine,
School of Mathematical Sciences,
University of Nottingham,
Nottingham NG7 2RD,
England,
Tomas Alarcon† and Philip K. Maini,
Centre for Mathematical Biology,
Mathematical Institute,
24-29 St Giles’,
Oxford OX1 3LB,
England.
ABSTRACT
Cancer is a complex disease in which a variety of phenomena
interact over a wide range of spatial and temporal scales. In
this article a theoretical framework will be introduced that is
capable of linking together such processes to produce a detailed
model of vascular tumour growth. The model is formulated
as a hybrid cellular automaton and contains submodels
that describe subcellular, cellular and tissue level features.
Model simulations will be presented to illustrate the effect
that coupling between these different elements has on the tumour’s
evolution and its response to chemotherapy.
1. INTRODUCTION
Cancer is a complex and insidious disease in which controls
designed to regulate growth and maintain homeostasis become
disrupted. It is frequently initiated by genetic mutations
that increase the net rate of cell division and lead to the formation
of a small avascular lesion. Successful angiogenesis (i.e.
the formation and ingrowth of a network of new blood vessels
to the tumour) is needed before the extensive and rapid growth
associated with vascularised tumours can occur. To confound
matters, the processes involved in tumour growth are interlinked
and act over a wide spectrum of spatial and temporal
scales: the spatial scales of interest range from the subcellular
level, to the cellular and macroscopic levels while the
timescales may vary from seconds (or less) for signal transduction
pathways to months for tumour doubling times.
The advent of increasingly sophisticated technology means
that it is now possible to collect experimental data associated
with the spatial and temporal scales of interest. This is creating
a demand for new theoretical models that have the capacity
to integrate such data in a meaningful manner and are
∗The authors gratefully acknowledge financial support provided by the
EPSRC under grants GR/5090067, GR/S72023/01 and AF/00067.
†The third author is currently based at the Bioinformatics Unit, Department
of Computer Science, University College London, Gower Street, London
WC1E 6BT, England.
able to address the fundamental problem of how phenomena
at different spatial scales are coupled.
In this paper we review our recent progress in developing
a mathematical model for studying vascular tumour growth
that is capable of integrating phenomena that act on different
scales [1]. Our theoretical framework extends earlier work
by Gatenby and coworkers [2] and links submodels which
describe processes operating on different spatial scales. In
section 2 we introduce our hybrid cellular automaton, before
presenting numerical results in section 3. The simulations illustrate
how the coupling between the submodels influences
the tumour’s evolution and its response to chemotherapy. We
conclude in section 4 with a summary of our results and a
discussion of possible directions for future research.
2. MODEL FORMULATION
In this section we introduce our multiscale model of vascular
tumour growth. It accounts for a variety of inter-related phenomena
that operate on vastly different space and time scales.
We consider a vasculature composed of a regular hexagonal
network embedded in a two-dimensional NxN lattice composed
of normal cells, cancer cells and empty space. Progress
through the cell cycle and the production of proteins such
as vascular endothelial growth factor (VEGF) that stimulate
angiogenesis are incorporated at the subcellular level using
ODE models. Cell-cell communication and competition for
resources are included at the cellular level through rules that
define our cellular automaton. At the tissue scale, reactiondiffusion
equations model the diffusion, production and uptake
of oxygen and VEGF: the vessels are regarded as sources
(sinks) of oxygen (VEGF) and the cells as sinks (sources)
of oxygen (VEGF). Blood flow and vascular adaptation are
also included at the tissue scale (see Figure 1). We impose
a pressure drop across the vasculature, assuming that blood
flows into the idealised “tissue” through a single inlet vessel
and drains out through a single outlet vessel. We use the
Poiseuille approximation and compute the flow rates through,
and pressure drops across, each vessel using Kirchoff’s laws.
The vessel radii are updated using a structural adaptation law
similar to that proposed by Pries et al. [3] (for details, see
[1]).
Thus themodel is formulated as a hybrid cellular automata,
with different submodels describing behaviour at the subcellular,
cellular and macroscopic (or vascular) levels (see Figure
1). Coupling between the different submodels is achieved in
several ways. For example, local oxygen levels which are determined
at the macroscale influence both progress through
the cell cycle and VEGF production at the subcellular level.
Conversely, the intracellular production of VEGF modulates
vascular adaptation at the macroscale and this, in turn, controls
oxygen delivery to the tissue. We stress that the submodels
we use simply illustrate how such a multiscale model can
be assembled: the framework we present is general, with considerable
scope for incorporatingmore realistic (i.e. complex)
submodels. This raises the important issue of how the level
of detail incorporated at each spatial scale influences the system’s
behaviour: this will form the basis of future research.
NUMERICAL RESULTS
3.1. Vascular adaptation influences tumour growth
In Figure 2 we present simulations that illustrate the importance
of accurately modelling blood flow through the tissue.
The upper panels correspond to a case for which the vessels
undergo structural adaptation and, hence, oxygen is distributed
nonuniformly across the tissue. The lower panels
show how the system evolves when oxygen is distributed uniformly
throughout the vessels (i.e. blood flow is identical in
all branches of the vasculature). We see that spatial heterogeneity
has a significant effect on the tumour’s dynamics and,that if the oxygen distribution is heterogeneous then the tumour
has “finger-like” protrusions similar to those observed
in invasive cancers. This structure arises here simply because
of the spatial heterogeneity in the nutrient distribution. Indeed,
closer inspection reveals that several parts of the tumour
have almost “broken away”. While this cannot actually happen
in the current model because cell motion is neglected, we
speculate that by allowing cell movement towards nutrientrich
regions, this may act as a mechanism for metastasis.
Fig. 2. Series of images showing the spatial distribution of
cells for growth in inhomogeneous (panel a), and homogeneous
environments (panels c). In panels (a) and (c) cancer
cells occupy white spaces and vessels occupy a hexagonal arrray
denoted by black spaces. The other black spaces denote
”empty spaces”. Panels (b) and (d) show the time evolution of
the number of (cancer) cells for the heterogeneous and homogeneous
cases, respectively: squares denote the total number
of cancer cells (proliferating + quiescent); diamonds denote
quiescent cells. Reproduced with permission from [1].
3.2. Impact of VEGF on the tumour’s growth dynamics
The simulation presented in Figures 3 to 5 shows how coupling
intracellular and macroscale phenomena can influence
the dynamics of both the vasculature and the tumour. In contrast
to the results depicted in Figure 2, where vessel adaptation
was independent of VEGF, in Figures 3 to 5 it is regulated
by local VEGF levels. Figures 3 and 4 show how the tumour’s
spatial composition evolves while figure 5 summarises its dynamics.
Since there is a single inlet (outlet) to the vasculature
located in the bottom left (top right) hand corner of the tissue,
the incoming blood flow and haematocrit become diluted as
they pass through the hexagonal lattice. This creates a het-
in this case, actually reduces the tumour burden. We note al
erogeneous oxygen distribution across the domain, with oxygen
levels being highest near the inlet and outlet. Over time,
the tumour cells proliferate and spread through the tissue towards
oxygen-rich regions. As they increase in number, their
demand for oxygen outstrips that available from the vasculature,
and quiescent regions form. These cells produce VEGF
which diffuses through the tissue (see Figures 3 and 4), stimulating
vessel adaptation and biasing blood flow towards low
oxygen regions. If the VEGF stimulus is weak then the vasculature
does not adapt quickly enough and the quiescent cells
die (this is what happens at early times in figure 5). VEGF
levels also decline and blood flow to the remaining tumour
cells rises, enabling them to increase in number until the demand
for oxygen once again exceeds that being supplied, and
so the cycle repeats, with pronounced oscillations in the number
of quiescent cells (see Figure 5). In order to highlight the
key role played by VEGF in creating these oscillations, also
presented in Figure 5 are the results of a simulation which
was identical in all respects except that vascular adaptation
was independent of VEGF (as per Figure 2). In both cases,
the tumours grow to similar sizes. However, when vascular
adaptation is independent of VEGF the evolution is monotonic,
the oscillations in the cell populations disappear and the
number of quiescent cells is consistently much lower. These
results show howcoupling between the different spatial scales
can effect not only the tumour’s growth dynamics but also the
proportion of proliferating and quiescent cells that it contains.Fig. 3. Series of plots showing how a small tumour introduced
into a vascular tissue at t = 0 has evolved at t = 30
(dimensionless time unit). While the oxygen and vessel profiles
remain unchanged from their initial configurations, the
tumour has increased in size and now contains quiescent cells
which produce trace amounts of VEGF. Reproduced with permission
from [4].Fig. 4. Series of plots showing how the simulation presented
in figure 3 has developed at t = 90. The tumour continues
to penetrate the tissue region. There are now enough quiescent
cells to elicit an angiogenic response. As a result, the
vasculature has been remodelled, with blood flow and oxygen
supply (haematocrit) being directed primarily towards the tumour
mass. Reproduced with permission from [4].
Response to chemotherapy
We now investigate how the system’s dynamics change when
a chemotherapeutic drug is introduced. We assume that the
drug is continuously administered to the vessels and, hence,
that its concentration at the inlet vessel is constant. We calculate
the drug concentration within and outside the vessels inFig. 5. Series of curves showing how, for the simulation in
figures 3 and 4 the numbers of proliferating (upper panel),
quiescent tumour cells (middle panel) and total number of tumour
cells (lower panel) change over time. While the number
of proliferating cells increases steadily, the number of quiescent
cells undergoes oscillations of increasing amplitude until
t ≈ 120. Thereafter, the tumour is sufficiently large that
the quiescent cells are never eliminated: quiescent cells that
die are replaced by proliferating cells that become quiescent.
The dot-dashed lines show the evolution of a tumour which is
identical except that its vasculature is not regulated by VEGF.
While both tumours reach similar equilibrium sizes, when
vascular adaptation is independent of VEGF the oscillations
in the cell populations disappear and the number of quiescent
cells is much lower. Reproduced with permission from [4].
4. DISCUSSION
We have presented a hybrid cellular automaton model of vascular
tumour growth and shown how it may be used to study
the manner in which interactions between subcellular, cellular
and macroscale phenomena affect the tumour’s growth dynamics
and its response to chemotherapy. We stress that the
submodels we have used to describe the different processes
are highly idealised and chosen simply to illustrate the potential
value of such a multiscale model as a predictive tool to
test experimental hypotheses and to integrate different types
of experimental data. There is considerable scope for incorporatingmore
realistic submodels and specialising the system
to describe specific tumour types. For example, we are currently
engaged in a large interdisciplinary project which aims
to build a virtual model of the early stages of colorectal cancer
(details at: http://www.integrativebiology.ox.ac.uk ).
Key challenges raised by our simulations that lie at the
heart of such integrative modelling concern the level of detail
incorporated at each spatial scale, the mathematical approaches
used and model validation. For example, in this
Series of curves showing how the tumour’s overall
dynamics change when it is treated with two blood-borne
chemotherapeutic agents that differ only in their extravasation
rates. For each simulation, we plot the numbers of proliferating
and quiescent tumour cells and the total number of tumour
cells evolve over time. Key: hθ = 0 (control, drug-free case,
as per figure 5), solid line; hθ = 90 (moderate drug), dashed
line; hθ = 100 (highly permeable and effective drug), dotted
line. Reproduced with permission from [4].
article we chose to use a combination of differential equations
and cellular automata to construct our virtual tumour.
It remains an open question whether the predicted behaviour
would change if we replaced our (subcellular) ODE models
with Boolean networks and/or the cellular automata with
agent-based models.
REFERENCES
[1] T. Alarc´on, H.M. Byrne, and P.K. Maini, “A muliple scale
model for tumour growth,” SIAM J. Multiscale Mod. &
Sim., vol. 3, pp. 440–475, 2005.
[2] A.A. Patel, E.T. Gawlinsky, S.K. Lemieux, and R.A.
Gatenby, “Cellular automaton model of early tumour
growth and invasion: the effects of native tissue vascularity
and increased anaerobic tumour metabolism,” J.
Theor. Biol., vol. 213, pp. 315–331, 2001.
[3] A.R. Pries, T.W. Secomb, and P. Gaehtgens, “Structural
adaptation and stability of microvascular networks: theory
and simulations,” Am. J. Physiol., vol. 275, pp.H349–
H360, 1998.
[4] H.M. Byrne, M.R. Owen, T. Alarc´on, J. Murphy, and
P.K. Maini (in press), “Modelling the response of vascular
tumours to chemotherapy: a multiscale approach,”
Math. Mod. Meth. Appl. Sci., 2005.
Cancer is a complex disease in which a variety of phenomena
interact over a wide range of spatial and temporal scales. In
this article a theoretical framework will be introduced that is
capable of linking together such processes to produce a detailed
model of vascular tumour growth. The model is formulated
as a hybrid cellular automaton and contains submodels
that describe subcellular, cellular and tissue level features.
Model simulations will be presented to illustrate the effect
that coupling between these different elements has on the tumour’s
evolution and its response to chemotherapy.
1. INTRODUCTION
Cancer is a complex and insidious disease in which controls
disrupted. It is frequently initiated by genetic mutations
that increase the net rate of cell division and lead to the formation
of a small avascular lesion. Successful angiogenesis (i.e.
the formation and ingrowth of a network of new blood vessels
to the tumour) is needed before the extensive and rapid growth
associated with vascularised tumours can occur. To confound
matters, the processes involved in tumour growth are interlinked
and act over a wide spectrum of spatial and temporal
scales: the spatial scales of interest range from the subcellular
level, to the cellular and macroscopic levels while the
timescales may vary from seconds (or less) for signal transduction
pathways to months for tumour doubling times.
The advent of increasingly sophisticated technology means
that it is now possible to collect experimental data associated
with the spatial and temporal scales of interest. This is creating
a demand for new theoretical models that have the capacity
to integrate such data in a meaningful manner and are
∗The authors gratefully acknowledge financial support provided by the
EPSRC under grants GR/5090067, GR/S72023/01 and AF/00067.
†The third author is currently based at the Bioinformatics Unit, Department
of Computer Science, University College London, Gower Street, London
WC1E 6BT, England.
able to address the fundamental problem of how phenomena
at different spatial scales are coupled.
In this paper we review our recent progress in developing
a mathematical model for studying vascular tumour growth
that is capable of integrating phenomena that act on different
scales [1]. Our theoretical framework extends earlier work
by Gatenby and coworkers [2] and links submodels which
describe processes operating on different spatial scales. In
section 2 we introduce our hybrid cellular automaton, before
presenting numerical results in section 3. The simulations illustrate
how the coupling between the submodels influences
the tumour’s evolution and its response to chemotherapy. We
conclude in section 4 with a summary of our results and a
discussion of possible directions for future research.
2. MODEL FORMULATION
In this section we introduce our multiscale model of vascular
tumour growth. It accounts for a variety of inter-related phenomena
that operate on vastly different space and time scales.
We consider a vasculature composed of a regular hexagonal
network embedded in a two-dimensional NxN lattice composed
of normal cells, cancer cells and empty space. Progress
through the cell cycle and the production of proteins such
as vascular endothelial growth factor (VEGF) that stimulate
angiogenesis are incorporated at the subcellular level using
ODE models. Cell-cell communication and competition for
resources are included at the cellular level through rules that
define our cellular automaton. At the tissue scale, reactiondiffusion
equations model the diffusion, production and uptake
of oxygen and VEGF: the vessels are regarded as sources
(sinks) of oxygen (VEGF) and the cells as sinks (sources)
of oxygen (VEGF). Blood flow and vascular adaptation are
also included at the tissue scale (see Figure 1). We impose
a pressure drop across the vasculature, assuming that blood
flows into the idealised “tissue” through a single inlet vessel
and drains out through a single outlet vessel. We use the
Poiseuille approximation and compute the flow rates through,
and pressure drops across, each vessel using Kirchoff’s laws.
The vessel radii are updated using a structural adaptation law
similar to that proposed by Pries et al. [3] (for details, see
[1]).
Thus themodel is formulated as a hybrid cellular automata,
with different submodels describing behaviour at the subcellular,
cellular and macroscopic (or vascular) levels (see Figure
1). Coupling between the different submodels is achieved in
several ways. For example, local oxygen levels which are determined
at the macroscale influence both progress through
the cell cycle and VEGF production at the subcellular level.
Conversely, the intracellular production of VEGF modulates
vascular adaptation at the macroscale and this, in turn, controls
oxygen delivery to the tissue. We stress that the submodels
we use simply illustrate how such a multiscale model can
be assembled: the framework we present is general, with considerable
scope for incorporatingmore realistic (i.e. complex)
submodels. This raises the important issue of how the level
of detail incorporated at each spatial scale influences the system’s
behaviour: this will form the basis of future research.
NUMERICAL RESULTS
3.1. Vascular adaptation influences tumour growth
In Figure 2 we present simulations that illustrate the importance
of accurately modelling blood flow through the tissue.
The upper panels correspond to a case for which the vessels
undergo structural adaptation and, hence, oxygen is distributed
nonuniformly across the tissue. The lower panels
show how the system evolves when oxygen is distributed uniformly
throughout the vessels (i.e. blood flow is identical in
all branches of the vasculature). We see that spatial heterogeneity
has a significant effect on the tumour’s dynamics and,that if the oxygen distribution is heterogeneous then the tumour
has “finger-like” protrusions similar to those observed
in invasive cancers. This structure arises here simply because
of the spatial heterogeneity in the nutrient distribution. Indeed,
closer inspection reveals that several parts of the tumour
have almost “broken away”. While this cannot actually happen
in the current model because cell motion is neglected, we
speculate that by allowing cell movement towards nutrientrich
regions, this may act as a mechanism for metastasis.
Fig. 2. Series of images showing the spatial distribution of
cells for growth in inhomogeneous (panel a), and homogeneous
environments (panels c). In panels (a) and (c) cancer
cells occupy white spaces and vessels occupy a hexagonal arrray
denoted by black spaces. The other black spaces denote
”empty spaces”. Panels (b) and (d) show the time evolution of
the number of (cancer) cells for the heterogeneous and homogeneous
cases, respectively: squares denote the total number
of cancer cells (proliferating + quiescent); diamonds denote
quiescent cells. Reproduced with permission from [1].
3.2. Impact of VEGF on the tumour’s growth dynamics
The simulation presented in Figures 3 to 5 shows how coupling
intracellular and macroscale phenomena can influence
the dynamics of both the vasculature and the tumour. In contrast
to the results depicted in Figure 2, where vessel adaptation
was independent of VEGF, in Figures 3 to 5 it is regulated
by local VEGF levels. Figures 3 and 4 show how the tumour’s
spatial composition evolves while figure 5 summarises its dynamics.
Since there is a single inlet (outlet) to the vasculature
located in the bottom left (top right) hand corner of the tissue,
the incoming blood flow and haematocrit become diluted as
they pass through the hexagonal lattice. This creates a het-
in this case, actually reduces the tumour burden. We note al
erogeneous oxygen distribution across the domain, with oxygen
levels being highest near the inlet and outlet. Over time,
the tumour cells proliferate and spread through the tissue towards
oxygen-rich regions. As they increase in number, their
demand for oxygen outstrips that available from the vasculature,
and quiescent regions form. These cells produce VEGF
which diffuses through the tissue (see Figures 3 and 4), stimulating
vessel adaptation and biasing blood flow towards low
oxygen regions. If the VEGF stimulus is weak then the vasculature
does not adapt quickly enough and the quiescent cells
die (this is what happens at early times in figure 5). VEGF
levels also decline and blood flow to the remaining tumour
cells rises, enabling them to increase in number until the demand
for oxygen once again exceeds that being supplied, and
so the cycle repeats, with pronounced oscillations in the number
of quiescent cells (see Figure 5). In order to highlight the
key role played by VEGF in creating these oscillations, also
presented in Figure 5 are the results of a simulation which
was identical in all respects except that vascular adaptation
was independent of VEGF (as per Figure 2). In both cases,
the tumours grow to similar sizes. However, when vascular
adaptation is independent of VEGF the evolution is monotonic,
the oscillations in the cell populations disappear and the
number of quiescent cells is consistently much lower. These
results show howcoupling between the different spatial scales
can effect not only the tumour’s growth dynamics but also the
proportion of proliferating and quiescent cells that it contains.Fig. 3. Series of plots showing how a small tumour introduced
into a vascular tissue at t = 0 has evolved at t = 30
(dimensionless time unit). While the oxygen and vessel profiles
remain unchanged from their initial configurations, the
tumour has increased in size and now contains quiescent cells
which produce trace amounts of VEGF. Reproduced with permission
from [4].Fig. 4. Series of plots showing how the simulation presented
in figure 3 has developed at t = 90. The tumour continues
to penetrate the tissue region. There are now enough quiescent
cells to elicit an angiogenic response. As a result, the
vasculature has been remodelled, with blood flow and oxygen
supply (haematocrit) being directed primarily towards the tumour
mass. Reproduced with permission from [4].
Response to chemotherapy
We now investigate how the system’s dynamics change when
a chemotherapeutic drug is introduced. We assume that the
drug is continuously administered to the vessels and, hence,
that its concentration at the inlet vessel is constant. We calculate
the drug concentration within and outside the vessels inFig. 5. Series of curves showing how, for the simulation in
figures 3 and 4 the numbers of proliferating (upper panel),
quiescent tumour cells (middle panel) and total number of tumour
cells (lower panel) change over time. While the number
of proliferating cells increases steadily, the number of quiescent
cells undergoes oscillations of increasing amplitude until
t ≈ 120. Thereafter, the tumour is sufficiently large that
the quiescent cells are never eliminated: quiescent cells that
die are replaced by proliferating cells that become quiescent.
The dot-dashed lines show the evolution of a tumour which is
identical except that its vasculature is not regulated by VEGF.
While both tumours reach similar equilibrium sizes, when
vascular adaptation is independent of VEGF the oscillations
in the cell populations disappear and the number of quiescent
cells is much lower. Reproduced with permission from [4].
4. DISCUSSION
We have presented a hybrid cellular automaton model of vascular
tumour growth and shown how it may be used to study
the manner in which interactions between subcellular, cellular
and macroscale phenomena affect the tumour’s growth dynamics
and its response to chemotherapy. We stress that the
submodels we have used to describe the different processes
are highly idealised and chosen simply to illustrate the potential
value of such a multiscale model as a predictive tool to
test experimental hypotheses and to integrate different types
of experimental data. There is considerable scope for incorporatingmore
realistic submodels and specialising the system
to describe specific tumour types. For example, we are currently
engaged in a large interdisciplinary project which aims
to build a virtual model of the early stages of colorectal cancer
(details at: http://www.integrativebiology.ox.ac.uk ).
Key challenges raised by our simulations that lie at the
heart of such integrative modelling concern the level of detail
incorporated at each spatial scale, the mathematical approaches
used and model validation. For example, in this
Series of curves showing how the tumour’s overall
dynamics change when it is treated with two blood-borne
chemotherapeutic agents that differ only in their extravasation
rates. For each simulation, we plot the numbers of proliferating
and quiescent tumour cells and the total number of tumour
cells evolve over time. Key: hθ = 0 (control, drug-free case,
as per figure 5), solid line; hθ = 90 (moderate drug), dashed
line; hθ = 100 (highly permeable and effective drug), dotted
line. Reproduced with permission from [4].
article we chose to use a combination of differential equations
and cellular automata to construct our virtual tumour.
It remains an open question whether the predicted behaviour
would change if we replaced our (subcellular) ODE models
with Boolean networks and/or the cellular automata with
agent-based models.
REFERENCES
[1] T. Alarc´on, H.M. Byrne, and P.K. Maini, “A muliple scale
model for tumour growth,” SIAM J. Multiscale Mod. &
Sim., vol. 3, pp. 440–475, 2005.
[2] A.A. Patel, E.T. Gawlinsky, S.K. Lemieux, and R.A.
Gatenby, “Cellular automaton model of early tumour
growth and invasion: the effects of native tissue vascularity
and increased anaerobic tumour metabolism,” J.
Theor. Biol., vol. 213, pp. 315–331, 2001.
[3] A.R. Pries, T.W. Secomb, and P. Gaehtgens, “Structural
adaptation and stability of microvascular networks: theory
and simulations,” Am. J. Physiol., vol. 275, pp.H349–
H360, 1998.
[4] H.M. Byrne, M.R. Owen, T. Alarc´on, J. Murphy, and
P.K. Maini (in press), “Modelling the response of vascular
tumours to chemotherapy: a multiscale approach,”
Math. Mod. Meth. Appl. Sci., 2005.
Tag :
HEALTH CARE
0 Komentar untuk "CANCER DISEASE"