Disease processes as hybrid dynamical systems
We investigate the use of hybrid techniques in complex processes of infectious diseases. Since predictive disease models in biomedicine require a multiscale approach for understanding the molecule-cell-tissue-organ-body interactions, heterogeneous methodologies are often employed for describing the different biological scales. Hybrid models provide effective means for complex disease modelling where the action and dosage of a drug or a therapy could be meaningfully investigated: the infection dynamics can be classically described in a continuous fashion, while the scheduling of multiple treatment discretely. We define an algebraic language for specifying general disease processes and multiple treatments, from which a semantics in terms of hybrid dynamical system can be derived. Then, the application of control-theoretic tools is proposed in order to compute the optimal scheduling of multiple therapies. The potentialities of our approach are shown in the case study of the SIR epidemic model and we discuss its applicability on osteomyelitis, a bacterial infection affecting the bone remodelling system in a specific and multiscale manner. We report that formal languages are helpful in giving a general homogeneous formulation for the different scales involved in a multiscale disease process; and that the combination of hybrid modelling and control theory provides solid grounds for computational medicine.
eptcsdata.tex
\providecommand{\firstpage}{152}
\providecommand{\publicationstatus}{Ezio Bartocci and Luca Bortolussi (Eds.): HSB 2012\\
EPTCS 92, 2012, pp.\ \pageref*{FirstPage}--\pageref{LastPage},
doi:\href{http://dx.doi.org/10.4204/EPTCS.92.11}{10.4204/EPTCS.92.11}}
\pdfoutput=1

hsb_camera_ready_ver2.tex
\documentclass[submission,copyright,creativecommons]{eptcs}
\providecommand{\event}{HSB 2012} % Name of the event you are submitting to
\usepackage{color}
\usepackage{array}
\usepackage{amsmath,amssymb}
\usepackage{graphicx}
%\usepackage{algorithm}
\usepackage{faktor}
%\usepackage{algorithmic}
\usepackage[caption=false,font=footnotesize]{subfig}

\newtheorem{thm}{Theorem}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
%\theoremstyle{definition}
\newtheorem{defn}[thm]{\ \\Definition}
%\theoremstyle{remark}
\newtheorem{rem}[thm]{\ \\Remark}

\newtheorem{example}[thm]{\ \\Example}
\renewcommand{\theequation}{\thesection.\arabic{equation}}

\title{Disease processes as hybrid dynamical systems}
\author{Pietro Li\`{o}
\institute{Computer Laboratory\\University of Cambridge\\Cambridge, UK}
\email{pl219@cam.ac.uk}
\and
Emanuela Merelli \qquad\qquad Nicola Paoletti
\institute{School of Science and Technology\\Computer Science Division, University of Camerino\\Camerino, IT}
\email{\quad emanuela.merelli@unicam.it \quad\qquad nicola.paoletti@unicam.it}
}
\def\titlerunning{Disease processes as hybrid dynamical systems}
\def\authorrunning{P. Li\`{o}, E. Merelli \& N. Paoletti}
\begin{document}
\maketitle

\begin{abstract}
We investigate the use of hybrid techniques in complex processes of infectious diseases. Since predictive disease models in biomedicine require a multiscale approach for understanding the molecule-cell-tissue-organ-body interactions, heterogeneous methodologies are often employed for describing the different biological scales. Hybrid models provide effective means for complex disease modelling where the action and dosage of a drug or a therapy could be meaningfully investigated: the infection dynamics can be classically described in a continuous fashion, while the scheduling of multiple treatment discretely. We define an algebraic language for specifying general disease processes and multiple treatments, from which a semantics in terms of hybrid dynamical system can be derived. Then, the application of control-theoretic tools is proposed in order to compute the optimal scheduling of multiple therapies. The potentialities of our approach are shown in the case study of the SIR epidemic model and we discuss its applicability on osteomyelitis, a bacterial infection affecting the bone remodelling system in a specific and multiscale manner. We report that formal languages are helpful in giving a general homogeneous formulation for the different scales involved in a multiscale disease process; and that the combination of hybrid modelling and control theory provides solid grounds for computational medicine.
\end{abstract}

\section{How many scales does it take to describe an epidemics?}
In this work we investigate multi-methodology approaches for the modelling of complex infectious diseases and for finding the best strategy of intervention, able to avoid or limit the spread of the disease and its effect on the affected organisms. Infectious diseases can be broadly classified into three groups (Fig.~\ref{fig:classification}): acute, latent persistent and chronic persistent. Acute diseases like the common cold, the Rhinovirus, the Yellow Fever, the Influenza or some strains of the Staphylococcus Aureus (Osteomyelitis) are characterized by a single disease episode after which they do not occur anymore. Latent persistent ones like the Herpes simplex, the Varicella-zoster or the Measles-SSPE arises also after the first disease episode and a non-infectious and latent period. Chronic persistent diseases (e.g.\ Hepatitis B, HIV, HTLV-1 leukemia, chronic Osteomyelitis) can protract their effects on the host organism for several years. 

Epidemic modelling is one of the most established tools for predicting the progress and the spread of a disease in large populations, which are commonly divided into compartments:%(Fig.~\ref{fig:classification} (b)):
\textbf{S}usceptible, \textbf{E}xposed, \textbf{I}nfected and \textbf{R}ecovered. According to the possible flows among such compartments, different variants arises, like the SIR (Susceptible$\rightarrow$Infected$\rightarrow$Recovered) model~\cite{kermack1932contributions}, SIS, SIRS, SEIR and SEIRS. Further compartments can be taken into account like the population of immune infants, usually indicated with M. Similar compartments can be identified not only at the human population level, but also at the cellular population level, where the pathogen typically acts by infecting susceptible cells. As a matter of fact, infections are characterized by multiscale dynamics that affect the organism at multiple levels in the biological hierarchy, as shown on the left side of Fig.~\ref{fig:multiple_disease}: at the \textit{intracellular} level, in case of infection of a cell by a pathogen; at the \textit{intercellular/cellular population} level, in case of infection of susceptible cells by infected cells; at the \textit{tissue} level, when multi-cellular ensembles are involved; at the \textit{organ/individual} level, when infection spreads to other parts of the body; and at the \textit{human population} level, if the disease is transmitted among individuals. A further scale above the human population one could be the consciousness level for which several psychological and cognitive behaviour models exists. Just consider that psychological states like fear or stress can even affect (negatively) the immune system: a demonstration that all the scales in the biological hierarchy are intimately connected to each other. Their inherent multiscale dynamics make the modelling of infective processes a challenging and intriguing field of research. Indeed, pathogens and infective agents generally act at the cellular scale, but they need to exploit mechanisms at the human population level to spread and survive in other host organisms.
%Recent evidences~\cite{dingli2006cancer} suggest a connection between cancer and infectious diseases according to which cancer cells behave like infectious agents and move from one host to another, in the case of canine transmissible venereal tumour. Generally cancer is not transmitted among individuals, although an important percentage of cancer deaths are due to infectious agents that, instead, enables the transmission at the human population level. In addition, even the non-infectious cancer is characterized by multiscale dynamics similar to infections, since metastases can occur and the disease spreads to different parts of the body.
%\begin{figure}
%\centering
%\subfloat[Classification of infectious diseases]{\includegraphics[width=\textwidth]{img/disease_class_bw.pdf}}\\
%\subfloat[Compartmental epidemic models]{\includegraphics[width=0.25\textwidth]{img/sir.pdf}}
%\caption{(a) shows a simplified classification of infectious diseases. The disease intensity follows a Gaussian curve approximately partitioned in three stages. In the first one, the infection starts but any symptom is reported. In the second one called ``disease episode'' symptoms are reported and the disease can be transmitted to other individuals. In the last stage, the disease is over if it is acute or latent persistent. However some diseases can become chronic.
%(b) summarizes the flows in different compartmental epidemic models. The compartment considered are \textit{S}  (Susceptible), \textit{E} (Exposed), \textit{I} (Infected) and \textit{R} (Recovered). Some of the most common models are SIS (S$\rightarrow$I$\rightarrow$S), SIR, SIRS, SEIR and SEIRS.}
%\label{fig:classification}
%\end{figure}
\begin{figure}
\centering
\includegraphics[width=\textwidth]{img/disease_class_bw.pdf}
\caption{Simplified classification of infectious diseases. The disease intensity follows a Gaussian curve approximately partitioned in three stages. In the first one, the infection starts but any symptom is reported. In the second one called ``disease episode'' symptoms are reported and the disease can be transmitted to other individuals. In the last stage, the disease is over if it is acute or latent persistent. However some diseases can become chronic.}
\label{fig:classification}
\end{figure}

\subsection{Formal languages, hybrid modelling and control in multiscale infectious diseases}
\begin{figure}
\centering
\includegraphics[width=0.7\textwidth]{img/multi_level_infection_method.pdf}
\caption{Modelling of multiscale infective dynamics. Infections involve the intracellular level (a cell infected by a pathogen); the intercellular level (infection among cells); the tissue level (infection among collections of cells); the organ level (infection among different parts of the body); and the human population level (transmission among individuals). Each biological scale can be implemented with different modelling techniques like Gillespie's SSA, Continuous Time Markov Chains (CTMC), ordinary/partial differential equations (ODE/PDE), agent-based models (ABM), or finite state machines (FSM). Such heterogeneous semantics can be wrapped by a common formal language.}
\label{fig:multiple_disease}
\end{figure}
%In past few decades, a large amount of literature has been produced so far in the field of infectious disease modelling \hl{citation needed!}. In particular, modelling epidemics has been one of the most valuable tools in answering to how disease spread in population, as well as providing useful suggestions to policy-makers to prevent and estimate the risk and control the outbreak of diseases. The development of computational methodologies and tools has provided major insights and explanation into the molecular and population processes that drive epidemics. 
We believe that hybrid modelling could help in unravelling the complexity of the multiscale dynamics occurring in such diseases. We refer to 'hybrid modelling' not just with its classical meaning, i.e.\ the modelling of those systems characterized by the co-existence of continuous and discrete dynamics, but also with a methodological meaning. In multiscale systems different scales are typically approached with different methodologies, thus leading to the co-existence of heterogeneous modelling techniques, in other words to a methodologically hybrid modelling approach. The cellular level is typically represented with agent-based or ordinary differential equation models, while the tissue and organ levels are often described using image-based finite element modelling (partial differential equations). In the context of disease modelling the determination of suitable intervention strategies, including drug and therapy administration, adds another level of description that affects multiple biological scales. 

This kind of integrative models are composed of single-scale models, describing the biological process at different characteristic space-time scales, and scale bridging models, which define how the single-scale models are coupled to each other~\cite{viceconti2011multiscale}. Higher level languages and formalisms can help in giving a general homogeneous formulation for the different scales in a multiscale biological system. In general, that formal description cannot be directly executed or simulated, but it is able to support multiple semantics (e.g. transition systems, differential equations and Markov chains). In particular process-algebraic languages are a formal notation initially developed for modelling software systems, but in the last decade it has been extensively used and extended in order to describe biological systems. In a seminal paper~\cite{cardelli2008processes}, Cardelli showed that a subset of stochastic CCS is powerful enough to encode both systems of reactions in the stoichiometric form and systems of ordinary differential equations. Additionally, a recent work of the authors~\cite{paoletti2011,	paoletti2011spec} demonstrates the effectiveness of a hybrid approach where a process-algebraic specification level is translated into a runnable stochastic agent-based model in the bone remodelling case study, while in~\cite{bartocci2012tcsb} the same case study is approached with different semantics according to the biological property to analyse.

Consequently, formal languages could represent a wrapping language able to homogeneously describe the different scales of a complex multiscale system, where each level can be instantiated into a runnable model according to the most suitable semantics. Figure~\ref{fig:multiple_disease} sketches the idea of the heterogeneous semantics associated to a multiscale model of infection through a common formal language. Furthermore, the importance of being hybrid in the classical sense is demonstrated by several biological evidences. For instance, genetic regulatory networks naturally exhibits hybrid dynamics, for which the continuous concentration of proteins is interrupted by the discrete switches dictated by changes in gene expression~\cite{batt2005qualitative,grosu2011cardiac}. Moreover in population biology it happens that some species are present in high concentration, so they can be approximately modelled by continuous variables; conversely, small populations are opportunely modelled as discrete stochastic variables~\cite{lapin2011shave}.

Another crucial aspect is related to the self-regulation and control of biological systems: in normal conditions, biological entities and functions are self-regulated and several multiscale control mechanisms naturally exists. One of the most striking example is the control operated by the immune system that protect the organism in case of disease. In fact, lymphocytes are able to detect the presence of a pathogen and produce an appropriate immune response by secreting immunoglobulins. However, in many cases of severe diseases the immune system cannot apply an effective control anymore. Drugs and therapies represent a form of control which is external to the organism and is extremely significant when dealing with models of diseases that can be limited by appropriate medical interventions. The control operated by medical therapies can be much more effective but needs to formulate control laws that take into account their (negative) impact on the organism due to possible side effects. While the drug administration strategy is determined at the human population level, the therapy affects multiple lower scales (in most of the cases, the cellular and subcellular scale). Other examples of multiscale control are the mechanical control at the organ and tissue level that regulates the functioning of bone cells in the bone remodelling process; or in an infection scenario, the cognitive behaviour at the human level that avoids or limits the interactions with other individuals possibly affected by the disease. A work by Bagnoli et al.~\cite{bagnoli2007risk} embeds this aspect of human behaviour into a SIR epidemic model by allowing individuals to perceive the risk of being infected by ill neighbours. Therefore, the use of control-theoretic tools seems promising for describing the biological mechanisms of multiscale self-regulation and self-adaptiveness, as well as for implementing externally imposed control mechanisms. %Effective control can be achieved on biological systems that are robust under certain conditions, and sensitive to others. According to Kitano~\cite{kitano2004biological}, \textit{``robustness is a property that allows a system to maintain its functions despite external and internal perturbations''}. Sensitivity analysis provides a measure for robustness since it studies the effects on model's output under perturbations from the expected input values. The design of (biological) control laws must take into account such properties.

%\iffalse
%Da mettere:
%\begin{itemize}
%\item \hl{controlli e scale}
%\item \hl{scala molecolare}
%\item \hl{controlli multipli e fitness}
%\item \hl{cosa succede al controllo con il disease e quindi}
%\item \hl{differenza tra controllo del sistema immunitario e terapie esterne}
%\item \hl{percezione del rischio (cognitive behaviour) come controllo, citare pietro}
%\end{itemize}
%
%We can distinguish between three types of control loops in a multiscale biological system:%, also depicted in Fig.~\ref{fig:multiscale_control}:
%\begin{itemize}
%\item \textit{intra-scale control}: entities at a specific scale perceive from each other information and signals useful to self-regulate their functions, e.g.\ cells are able to detect their neighbourhood, according to which drive their development;
%\item \textit{bottom-up control}: the dynamics at a specific scale is controlled by signals coming at lower scales, e.g.\ tissue homoeostasis can be seen as the objective function to be optimized by the controller at the lower scale (the cells);
%\item \textit{top-down control}: the dynamics at a specific scale is controlled by signals coming at higher scales, e.g.\ the bone remodelling effort at cellular level that is regulated by the mechanical stresses at tissue level, or the structure of the cardiovascular system that is controlled by the individual's metabolism and weight. In this case, higher scales are part of the environment of the considered scale.
%\end{itemize}
%\fi
%
%\begin{figure}
%\centering
%\includegraphics[width=0.5\textwidth]{img/multiscale_control_simple.pdf}
%\caption{Multiscale control in biological systems. The component to control is placed at scale $i$. According to the scale of the controller component three different control laws arises: \textit{intra-scale} (controller at scale $i$); \textit{bottom-up} (controller at scale $i - 1$); and \textit{top-down} (controller at scale $i + 1$).}
%\label{fig:multiscale_control}
%\end{figure}

In this work we present a general model for complex multiscale infectious diseases and for scheduling optimal medical treatments. We employ a hybrid modelling approach both in the classical and in the methodological sense. Indeed, we formulate the disease process by means of a high-level process algebra called D-CGF, from which a hybrid (in the classical sense) semantics can be derived. This semantics is given in terms of a hybrid dynamical system, where the cell population and the infection dynamics are described in a continuous fashion, while the dosage of multiple therapies is implemented as discrete switches. Then, we propose the application of model predictive control (MPC) tools for computing the optimal scheduling of multiple therapies. The potentialities of our approach are shown in the case study of the SIR epidemic model and we discuss its applicability on the case study of osteomyelitis, a bacterial infection affecting the bone remodelling system.

The paper is organized as follows. Section~\ref{sect:general_model} presents the algebraic language D-CGF for generic disease processes and its semantics in terms of hybrid dynamical systems, by means of the SIR example. Section~\ref{sect:sim} reports the optimal controlled solutions for scheduling multiple therapies in the SIR model and show how this approach can be useful to the osteomyelitis case study. Discussions and conclusions are given in Section~\ref{sect:conclusions}. 

\section{General process-algebraic formulation of infectious diseases}\label{sect:general_model}
In this part we propose a generic formulation of complex infectious diseases, by defining a variant of the Chemical Ground Form (CGF) stochastic process algebra. Stochastic process algebras extend classical ones with quantitative information in form of action rates and have associated a Continuous Time Markov Chains (CTMC) semantics from which a continuous and less computationally demanding approximation in terms of ODE can be derived. The multiplicity of semantics it supports has made formal process-algebraic languages one of the most used tools in biological modelling. A work by Cardelli~\cite{cardelli2008processes} shows how several encodings can be implemented among specifications in the CGF algebra (a subset of stochastic CCS), systems of reactions in the stoichiometric form and ODE systems. Bortolussi and Policriti~\cite{bortolussi2009hybrid} extend this work by defining a semantics in terms of hybrid automata~\cite{henzinger1996theory}, defining the requisites according to which a set of processes can be treated as discrete control states, while others have the usual continuous interpretation. Several other methods for associating hybrid semantics to quantitative process algebras have been proposed, among which~\cite{galpin2009hype,bortolussi2010hybrid}.

The formal modelling language employed in this paper is called \textit{Disease Chemical Ground Form (D-CGF)}, and it is a variant of CGF for describing complex disease processes. Differently from~\cite{norman2003developing}, where a stochastic process algebra and its continuous semantics are used to describe a classical epidemic model, we attempt to provide a framework for more generic diseases and, although the syntax of D-CGF does not radically differ from that of CGF, we define a novel constructive procedure for deriving a semantics in terms of hybrid dynamical systems, by distinguishing the following two kinds of processes:
\begin{itemize}
\item \textit{Individuals}: standard CGF processes that belong to some species and that collectively represent the population in the disease model. Different species can describe for instance the various compartments of an epidemic model (susceptible, infected, and etc.), a pathogen, or a particular mutation. Individuals are interpreted as continuous variables in the hybrid semantics.
\item \textit{Therapies}: processes for modelling interventions on the disease scenario. They could represent for instance the dosage of a drug or a change in the environment that \textit{discretely} alter the disease dynamics. Based on the hybrid semantics of CGF~\cite{bortolussi2009hybrid}, therapy processes are subject to some restrictions in order to be interpreted as the discrete switches of the hybrid dynamical system.
\end{itemize}
D-CFG will be illustrated by means of an epidemic example, even if this language is suitable to describe also models of virus infection, infectious cancers or other kinds of diseases where discrete intervention policies has to be modelled.
\subsection{The Disease Chemical Ground Form}
Here we present the \textit{Disease Chemical Ground Form (D-CGF)}, a variant of the Chemical Ground Form (CGF)~\cite{cardelli2008computational}) targeted to describe complex models of disease in a process-algebraic fashion. The syntax of \textit{D-CGF} is given by the grammar in Table~\ref{tbl:dcgfSyntax}. A D-CGF model is given by a set of species $S$, an initial population $P$, a set of therapies $T$, and an initial combination of therapies $C$. The first notable difference from the algebra CGF is that in D-CGF models two disjoint sets of processes are distinguished:
\begin{itemize}
\item the set of species (a set of individual species definitions), the individual species and the population (the parallel composition of individuals); and
\item the set of therapy definitions, the therapy and the combination of therapies (the parallel composition of therapies).
\end{itemize}
\begin{table}
\centering
\begin{small}
\begin{tabular}{rll}
$S::=$&$\mathbf{0} \ | \ X = I, S$ & Set of species definition\\
$I::=$&$\mathbf{0} \ | \ \pi . P + I$ & Individual species\\
$P::=$&$\mathbf{0} \ | ( X \| P)$ & Population\\
$\pi::=$&$\tau^r \ | \ ? x^r \ | \ ! x^r$ & Actions\\
&&\\
$T::=$&$\mathbf{0} \ | \ U = R, T$ & Set of therapy definition\\
$R::=$&$\mathbf{0} \ | \ \pi . C + R$ & Therapy\\
$C::=$&$\mathbf{0} \ | ( U \| C)$ & Combination of therapies\\
&&\\
$D-CGF::=$&$(S, P, T, C)$ & A \textit{D-CGF} model
\end{tabular}
\end{small}
\caption{Syntax of D-CGF.}\label{tbl:dcgfSyntax}
\end{table}
Individual species and therapies are defined by the alternative composition ($+$) of action-prefixed terms. Actions $\pi \in \Pi$ are indexed by a quantity $r$ corresponding to the rate of the exponential distribution that determines their duration, and are the usual ones: $\tau^r$ (internal action), $? x^r$ (input action) and $! x^r$ (output action). The motivation for separating these two sets of processes is related to the hybrid semantics of D-CGF: species are interpreted as continuous variables, while therapies as discrete control states. In addition we will need to put some restrictions on the therapy processes, for formally justifying this separation. In the remainder of this section, a simple epidemic model will serve as the running example for illustrating our approach.
\subsubsection{D-CGF by example}
In this part, we illustrate some of the technical features of the algebra by means of an epidemic example: the SIR (Susceptible$\rightarrow$Infected$\rightarrow$Recovered) model~\cite{kermack1932contributions}. It describes the dynamics at the epidemiological level of a population consisting of individuals susceptible to the disease ($S$), those infected ($I$) and those recovered ($R$). Here we consider an open-population (i.e.\ with births and deaths) SIR model, given by the list of reactions in Table~\ref{tbl:sirmodel}.
\begin{table}
\centering
\begin{small}
\begin{tabular}{rl}
$S \rightarrow_b S + S$&\\
$I \rightarrow_b I + S$&Birth of susceptible\\
$R \rightarrow_b R + S$&\\
\noalign{\smallskip}
$S \rightarrow_\mu \emptyset$&Death of susceptible\\
$I \rightarrow_\mu \emptyset$&Death of infected\\
$R \rightarrow_\mu \emptyset$&Death of recovered\\
\noalign{\smallskip}
$S + I \rightarrow_\beta I + I$&Infection of susceptible\\
\noalign{\smallskip}
$I \rightarrow_\nu R$&Infected becomes recovered
\end{tabular}
\end{small}
\caption{Reactions in the SIR epidemic model with open population and no disease-dependent death rate.}
\label{tbl:sirmodel}
\end{table}

According to the conversion rules for CGF, we can translate this list of reactions in a D-CGF model $\mathcal{M}=(S, P, \mathbf{0}, \mathbf{0})$ with no therapies.
{\small
\[
S:\ 
\begin{array}{rl}
S = & \tau_{S_1}^b . (S \| S) \ + \ \tau_{S_2}^\mu . \mathbf{0}\ + \ ? i^\beta . I\\
I = & \tau_{I_1}^b . (I \| S) \ + \ \tau_{I_2}^\mu . \mathbf{0}\ + \ ! i^\beta . I\ + \ \tau_{I_3}^\nu . R\\
R = & \tau_{R_1}^b . (R \| S) \ + \ \tau_{R_2}^\mu . \mathbf{0}
\end{array}
\qquad
P: \ P_0,
\]
}\noindent where $P_0$ is the initial population. Given an action $\pi \in \Pi$, we denote with $react(\pi)$ the multiset of species consumed by $\pi$, with $prod(\pi)$ the multiset of those produced by $\pi$. $\#(X,P)$ denotes the number of $X$ occurring in the population $P$. $\Delta (\pi, X) = \#(X,prod(\pi)) - \#(X,react(\pi))$ denotes the net variation of $X$ due to $\pi$. For instance in the above example, $react(i) = \lbrace S, I\rbrace$, $prod(i) = \lbrace I, I\rbrace$, $\Delta (i, S) = -1$ and $\Delta (i, I) = 1$.
%\iffalse
%\hl{We restrict D-CGF terms in order to consider birth, death and mutation actions, according to the following definition.}
%\begin{defn}[Allowed D-CGF actions]
%The set of allowed D-CGF actions contains the actions $\pi \in \Pi$ for which there exist individuals $X$ and $Y$ ($X \neq Y$) such that $\pi$ supports at least one of the following forms:
%\begin{itemize}
%\item \textit{Growth action for $X$:} $\Delta (\pi, X) > 0 \wedge react(\pi) \subseteq prod(\pi)$.
%\item \textit{Death action for $X$:} $\Delta (\pi, X) < 0 \wedge prod(\pi) \subseteq react(\pi)$. 
%\item \textit{Mutation action of $X$ into $Y$:} $\Delta (\pi, X) < 0 \wedge Y \in prod(\pi) \setminus react(\pi)$.
%\end{itemize}
%\end{defn}
%The growth formula requires that there is a net increase of $X$ and that the reactants of $\pi$ are a subset of its products, i.e.\ none of the reactants are consumed by $\pi$. The death formula tells that there is a net loss of $X$, but that the products of $\pi$ are a subset of its reactants, i.e.\ no new species are introduced by the reaction (no mutation). Finally, in a mutation action there is a net loss of $X$, but a new species $Y$ is introduced by $\pi$, indicating that $X$ has been mutated into $Y$. These basic rules are enough to model generic disease dynamics with births and deaths of individuals, and with mutations describing for instance a susceptible individual becoming infected. In the SIR model, $\tau_{S_1}$, $\tau_{I_1}$ and $\tau_{R_1}$ are growth actions for $S$; $\tau_{S_2}$, $\tau_{I_2}$ and $\tau_{R_2}$ are death actions for $S$, $I$ and $R$, respectively; $i$ is a mutation action of $S$ into $I$; and $\tau_{I_3}$ is a mutation of $I$ into $R$.
%\fi

The procedure for extracting a system of ODEs requires to build the so-called \textit{stoichiometric matrix} $M$, which has one row for each species and one column for each action. A cell $M[X,\pi]$ takes the net variation of $X$ due to $\pi$: $\Delta (\pi, X)$. Then a \textit{rate vector} $\phi$ has to be defined in the following way. Let $r(\pi)$ denote the rate of action $\pi$. Then for each action $\pi \in \Pi$,
{\small
\[
  \phi[\pi] = \left\{ 
  \begin{array}{l l}
    0 & \quad \text{if } react(\pi)=\emptyset\\
    r(\pi)X & \quad \text{if } react(\pi)=\lbrace X\rbrace\\
    r(\pi)XY & \quad \text{if } react(\pi)=\lbrace X, Y\rbrace\\
    r(\pi)X(X-1) & \quad \text{if } react(\pi)=\lbrace X, X\rbrace\\
  \end{array} \right.
\]
}\noindent Finally, the differential equations are given by $\mathbf{\dot{X}} = M \cdot \phi$. The stoichiometric matrix, the rate vector and the resulting ODEs for the SIR example are shown below.
{\footnotesize
\[
M = \bordermatrix{~ &\tau_{S_1} & \tau_{I_1} & \tau_{R_1} & \tau_{S_2} & \tau_{I_2} & \tau_{R_2} & \tau_{I_3} & i \cr
S& 1 & 1 & 1 & -1 & 0 & 0 & 0 & -1 \cr
I& 0 & 0 & 0 & 0 & -1 & 0 & -1 & 1 \cr
R& 0 & 0 & 0 & 0 & 0 & -1 & 1 & 0 \cr}
\qquad
\phi = \bordermatrix{~ &\cr
\tau_{S_1}&  bS\cr
\tau_{I_1}&  bI\cr
\tau_{R_1}&  bR\cr
\tau_{S_2}&  \mu S\cr
\tau_{I_2}&  \mu I\cr
\tau_{R_2}&  \mu R\cr
\tau_{I_3}&  \nu I\cr
i&  \beta SI \cr
}
\qquad
\begin{matrix}
\dot{S} = bN - \beta S I - \mu S\\
\dot{I} = \beta S I - \mu I - \nu I\\
\dot{R} = \nu I - \mu R\\
(N = S + I + R)
\end{matrix}
\]
}

Now we modify the SIR example in order to include two therapies $T1$ and $T2$ that makes immune the susceptible population and increase the mutation rate from infected to recovered, respectively. $T1_{on}$ ($T2_{on}$) and $T1_{off}$ ($T2_{on}$) denote the therapy being administered or not, respectively. Therefore, $T1$ can be thought as a vaccination that makes immune the susceptible individuals, while $T2$ as a therapy for combating the course of the disease. The therapy-specific reactions introduced are listed in Table~\ref{tbl:sirmodel1}.
\begin{table}
\centering
\begin{small}
\begin{tabular}{rl}
$S + T1_{on}\rightarrow_\rho R + T1_{on}$&Susceptible becomes recovered with therapy 1\\
\noalign{\smallskip}
$I + T2_{on}\rightarrow_{k} R + T2_{on}$&Infected becomes recovered with therapy 2\\
\noalign{\smallskip}
$T1_{off} \rightarrow_{r1_{on}} T1_{on}$& Therapy 1 switched on\\
$T1_{on} \rightarrow_{r1_{off}} T1_{off}$& Therapy 1 switched off\\
$T2_{off} \rightarrow_{r2_{on}} T2_{on}$& Therapy 2 switched on\\
$T2_{on} \rightarrow_{r2_{off}} T2_{off}$& Therapy 2 switched off\\
\end{tabular}
\end{small}
\caption{Additional therapy-specific reactions in the SIR epidemic model with open population and medical treatments at Table~\ref{tbl:sirmodel}.}
\label{tbl:sirmodel1}
\end{table}
In this case, the associated D-CGF model $\mathcal{M}=(S, P, T, C)$ has been extended with the therapy processes and is defined as follows:
{\small
\[
S: \ \begin{array}{rl}
S = & \tau_{S_1}^b . (S \| S) \ + \ \tau_{S_2}^\mu . \mathbf{0}\ + \ ? i^\beta . I \ + \ ? j^\rho . R\\
I = & \tau_{I_1}^b . (I \| S) \ + \ \tau_{I_2}^\mu . \mathbf{0}\ + \ ? \tau_{I_3}^\nu . R \ + \ ! i^\beta . I\ + \ ? h^{k} . R\\
R = & \tau_{R_1}^b . (R \| S) \ + \ \tau_{R_2}^\mu . \mathbf{0}
\end{array}
\qquad
P:\ P_0
\]\[
T: \ \begin{array}{rl}
T1_{off} = & \tau_{1on}^{r1_{on}} . T1_{on}\\
T1_{on} = & ! j^\rho . T1_{on} \ + \  \tau_{1off}^{r1_{off}} . T1_{off}\\
T2_{off} = & \tau_{2on}^{r2_{on}} . T2_{on}\\
T2_{on} = & ! h^{\nu + k} . T2_{on} \ + \ \tau_{2off}^{r2_{off}} . T2_{off}\\
\end{array}
\qquad
C: \ T1_{off} \ \| \ T2_{off}
\]
}
\subsection{Hybrid semantics of D-CGF}\label{subsect:hybridsemantics}
In the following we will give the definition of switching therapy (ST), and show that STs constitute the discrete switches in the hybrid dynamical system semantics where the different species are the continuous variables. The definition of switching therapy is broadly inspired by the notion of Control Automata~\cite{bortolussi2009hybrid}, i.e.\ the control structure that can be identified in the hybrid automata-based semantics of a CGF model. We also provide a constructing procedure for determining the STs associated to a D-CGF model and therefore its semantics in terms of hybrid dynamical systems.

%Before, some preliminary notions must be introduced. First, the structure of \textit{connection graph} that characterizes the links between reactants and products of D-CGF actions. 
%\begin{defn}[Connection graph]
%Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model. The connection graph $\mathcal{G}=(V,E)$ associated to $\mathcal{M}$ is a directed graph whose vertices are the types of species and therapies, $V = S \cup T$; and whose edges connect reactant processes to product processes, $E = \lbrace (X,Y) \subseteq V \times V \ | \ \exists \pi \in \Pi . X \in react(\pi) \wedge Y \in prod(\pi) \rbrace$.
%\end{defn}
%The connection graph for the SIR example with therapy is shown in Figure~\ref{fig:connectionGraph}. Given a set of process terms $S$, we denote with $\mathcal{G}[S]$ the sub-graph of $G$ restricted to the set of nodes $S$.
%\begin{figure}
%\centering
%\includegraphics[width=0.3\textwidth]{img/connection_graph.pdf}
%\caption{Connection graph for the SIR example. Edges connects reagents and products in a D-CGF action. The separation of species processes and therapy ones is for illustrative purposes and is not defined in the connection graph structure.}
%\label{fig:connectionGraph}
%\end{figure}

In order to interpret the set therapy terms $T$ as a set of discrete switches, we need to identify collections of terms in which exactly one term is active in every combination of therapies reachable from the initial one. We call such a collection a \textit{Switching Therapy (ST)}. The intuition is that a switching therapy models a discrete component whose terms represent its internal states, since exactly one of them must be active at each step. We give the definition of ST and of well-formed therapies for imposing additional restrictions on the construction of therapy terms in order to ensure a correct discrete interpretation of them.

\begin{defn}[Switching Therapy]\label{def:st}
Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model, and $T_i \subseteq T$ a set of therapy terms. $T_i$ is a \textit{Switching Therapy (ST)} if the following conditions hold:
\begin{enumerate}
%\item The connection graph $\mathcal{G}[T_i]$ is strongly connected.
\item Exactly one of the terms $U \in T_i$ is active in the initial combination of therapies: $\#(T_i, C) = \sum_{U \in T_i} \#(U, C) = 1$.
\item Each action must conserve the concentration of terms in $T_i$ and cannot involve more than one reagent in $T_i$: $\forall \pi \in \Pi. (\#(T_i, react(\pi))=\#(T_i, prod(\pi)) \leq 1)$.
%\item For each action $\pi \in \Pi$ such that $react[\pi]=\lbrace U_1, U_2\rbrace$, with $U_1 \in T_i$ and $U_2 \in T \setminus T_i$, $U_1 \not\in prod(\pi) \implies U_2 \in prod(\pi)$.
\item For each action $\pi \in \Pi$ such that $U_1 \in react(\pi)$ and $U_2 \in prod(\pi)$, with $U_1, U_2 \in T_i$, $U_1 \neq U_2 \implies (react(\pi) = \lbrace U_1\rbrace \wedge prod(\pi) = \lbrace U_2\rbrace)$.
\end{enumerate}
\end{defn}
Conditions 1+2 ensure that exactly one term of the ST is active in every reachable configuration, since condition 1 requires that the initial concentration of $T_i$-terms must be equal to one and condition 2 tells that such concentration is conserved. Additionally, condition 2 implies that there are no actions able to modify the concentration of $T_i$-terms and thus that the internal state of the ST $T_i$ can be changed only by actions having a $T_i$-term in its reactants. Finally, condition 3 requires that every action causing the switch of a $T_i$-term from $U_1$ to $U_2$ must be an internal action of $U_1$, of the form $U_1 = \ldots + \tau.U_2 + \ldots$.
%therapy terms belonging to different STs cannot switch simultaneously.
\begin{defn}[Well-formed therapies]\label{def:wellformed}
Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model and let $C=U_1 \| \ldots  \| U_n$, with $n \geq 1$. $T$ is a set of well-formed therapies if there exists a partition $\mathcal{T} = \lbrace T_1, \ldots , T_n \rbrace$ of $T$ such that $T_i \in \mathcal{T}$ is a switching therapy and $U_i \in T_i$, for all $i=1, \ldots, n$.
\end{defn}
We present a method for extracting the switching therapies and show how the semantics in terms hybrid dynamical systems can be constructed. Given a D-CGF model $\mathcal{M} = (S,P,T,C)$, the procedure exploits the stoichiometric matrix $M$, with which it is easy to check the conservation of switching terms and condition 3 of Definition~\ref{def:st}. Some necessary conditions has to be formulated on $M$ for ensuring the well-formedness of $T$.
\begin{defn}[Necessary conditions for well-formedness]\label{def:necessary}
Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model and $M$ be the associated stoichiometric matrix. Then, the following conditions are necessary for the well-formedness of $T$
\begin{enumerate}
\item The elements of $M$ restricted to $T$ take values in $\lbrace -1, 0, 1 \rbrace$: $\forall \pi \in \Pi. \ \forall U \in T. \ M_{|T}[U, \pi] \in \lbrace -1, 0, 1 \rbrace$.
\item Conservation of $T$-terms: $\forall \pi \in \Pi. \ \sum_{U \in T} M[U, \pi] \ = 0$;
\item Exclusive switch (1): $\forall \pi \in \Pi. \ \exists_{\leq 1} U \in T . \ M[U, \pi] = -1$. This condition and the conservation condition additionally imply that it exists exactly one $V \in T$ such that $M[V, \pi] = 1$.
\item Exclusive switch (2): $\forall \pi \in \Pi. \ (\exists U \in T . \ M[U, \pi] = -1 \implies (\forall X \in S . M[X, \pi]=0 \wedge \pi \text{ internal}))$. 
\end{enumerate}
\end{defn}

The stoichiometric matrix $M$ associated to the SIR example with therapies is given below.
{\footnotesize
\[
M = \bordermatrix{~ &\tau_{S_1} & \tau_{I_1} & \tau_{R_1} & \tau_{S_2} & \tau_{I_2} & \tau_{R_2} & \tau_{I_3} &\tau_{1on} & \tau_{1off} & \tau_{2on} & \tau_{2off} & i & j & h \cr
S& 1 & 1 & 1 & -1 & 0 & 0 &  0 & 0 & 0 & 0  &  0 & -1 & -1 & 0\cr
I& 0 & 0 & 0 & 0 & -1 & 0 &  -1 & 0 & 0 & 0  &  0 & 1 & 0 & -1\cr
R& 0 & 0 & 0 & 0 & 0 & -1 &  1 & 0 & 0 & 0  &  0 & 0 & 1 & 1\cr\cr
T1_{off}& 0 & 0 & 0 & 0 & 0 & 0 &  0 &-1 & 1 & 0  &  0 & 0 & 0 & 0\cr
T1_{on}& 0 & 0 & 0 & 0 & 0 & 0 &  0 & 1  & -1  & 0  &  0 & 0 & 0 & 0\cr
T2_{off}& 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0  & 0  & -1 &  1 & 0 & 0 & 0\cr
T2_{on}& 0 & 0 & 0 & 0 & 0 & 0 &  0 & 0  & 0  & 1  & -1 & 0 & 0 & 0\cr}
\]
}

It is easy to check that the conditions listed above are met in $M$.%, since for each $U \in T$, $M[U,\pi]$ takes -1, 0 or 1; for each action, the $T$-terms are conserved; for each action $\pi$, there is at most one couple of $T$-terms $(U_1,U_2)$ such that $M[U_1,\pi] = -1$ and $M_{|T}[U_2,\pi] = 1$, and for those actions $\pi$ it holds that $\pi$ is internal and $\forall X \in S . M[X, \pi]=0$.
Such conditions on the stoichiometric matrix are not sufficient to ensure the well-formedness of the set of therapies $T$. Indeed, Def.~\ref{def:necessary} corresponds to stating that the whole set $T$ meets conditions 2+3 in the definition of switching therapy, but it is not possible from the only stoichiometric matrix to prove also condition 1 (exactly one term in $T$ is active in the initial configuration). Moreover recalling that $T$ is well-formed if it can be partitioned into a set of switching therapies, such necessary conditions are defined on the entire set of therapies $T$, thus considering only the trivial partition $\lbrace T \rbrace$.

The algorithm for extracting the STs consists in building a graph $\mathcal{G} = (T, E)$ called \textit{ST-graph}, where $T$ is the set of process terms and arcs connect couples $(U_1,U_2)$ of therapy terms that are involved in a switch, i.e.\ $M[U_1, \pi] = -1$ and $M_{|T}[U_2, \pi] = 1$ for some $\pi$.

\begin{defn}[ST-graph]
Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model and $M$ be the associated stoichiometric matrix. The ST-graph $\mathcal{G_M}=(T,E)$ associated to $\mathcal{M}$ is a directed graph whose vertices are the $T$-terms of $\mathcal{M}$ and $E=\lbrace (U_1,U_2) \in T \times T \ | \ \exists \pi \in \Pi . \ (M[U_1, \pi] = -1 \wedge M_{|T}[U_2, \pi] = 1)\rbrace$.
\end{defn}

%Now it is possible to check the well-formedness of $T$ by considering the set of strongly connected components (SCC) of $\mathcal{G_M}$, $\mathcal{C(G_M)}$ (see Proposition~\ref{prop:wellform}). In particular, if the stoichiometric matrix $M$ meets the necessary conditions in Def.~\ref{def:necessary} and the set of SCCs $\mathcal{C(G_M)}=\lbrace G_1, \ldots, G_n\rbrace$ is such that its components $G_i$ are mutually disconnected and for each $G_i$ there is exactly one term active in the initial configuration, then $T$ is well-formed and the nodes of each $G_i$ form a switching therapy.
Now it is possible to check the well-formedness of $T$ by considering the set of connected components of $\mathcal{G_M}$, $\mathcal{C(G_M)}$ (see Proposition~\ref{prop:wellform}). In particular, if the stoichiometric matrix $M$ meets the necessary conditions in Def.~\ref{def:necessary} and the set of connected components $\mathcal{C(G_M)}=\lbrace G_1, \ldots, G_n\rbrace$ is such that for each $G_i$ there is exactly one term active in the initial configuration, then $T$ is well-formed and the nodes of each $G_i$ form a switching therapy.
%Note that the sets of nodes of the SCCs of an ST-graph form, by definition, a partition of $T$. This is necessary to guarantee that, according to Def.~\ref{def:wellformed}, a set well-formed therapies can be partitioned into a set of switching therapies. In addition, we have mentioned that the elements of $\mathcal{C(G_M)}$ must be mutually disconnected, i.e.\ given any two components $G_i, G_j \in \mathcal{C(G_M)}$ there cannot exist any arc in the ST-graph connecting them. This requirement ensures that Condition 2 of Def.~\ref{def:st} about the conservation of $T_i$-terms is met.
Note that the sets of nodes of the connected components in an ST-graph form, by definition, a partition of $T$. This is necessary to guarantee that, according to Def.~\ref{def:wellformed}, a set well-formed therapies can be partitioned into a set of switching therapies. In addition, the therapy terms in each $G_i$ meet the Condition 2 of Def.~\ref{def:st} about the conservation of terms in a switching therapy, because the connected components $G_i$ are (of course) mutually disconnected, i.e.\ there cannot exist any arc in the ST-graph connecting them that is, no terms can flow between them.

\begin{prop}\label{prop:wellform}
%Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model, $M$ be the associated stoichiometric matrix, $\mathcal{C(G_M)} = \lbrace G_1 = (T_1,E_1), \ldots , G_n  = (T_n,E_n) \rbrace$ be the set of strongly connected components of the associated ST-graph $\mathcal{G_M}$ and $\mathcal{T}=\lbrace T_1, \ldots, T_n \rbrace$. $T$ is well-formed and each $T_i \in \mathcal{T}$ is a switching therapy if the following conditions hold:
Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model, $M$ be the associated stoichiometric matrix, $\mathcal{C(G_M)} = \lbrace G_1 = (T_1,E_1), \ldots , G_n  = (T_n,E_n) \rbrace$ be the set of connected components of the associated ST-graph $\mathcal{G_M}$ and $\mathcal{T}=\lbrace T_1, \ldots, T_n \rbrace$. $T$ is well-formed and each $T_i \in \mathcal{T}$ is a switching therapy if the following conditions hold:
\begin{enumerate}
\item $M$ meets the necessary requirements at Def.~\ref{def:necessary};
%\item $G_1, \ldots G_n$ are mutually disconnected in $\mathcal{G_M}$;
\item for each $T_i \in \mathcal{T}$, exactly one of the terms in $T_i$ is active in the initial combination of therapies: $\#(T_i, C) = 1$; and
\item for each $T_i \in \mathcal{T}$ and for each action $\pi \in \Pi$, the number of $T_i$-terms in $react(\pi)$ is at most one: $\#(T_i, react(\pi)) \leq 1$.
\end{enumerate}
\end{prop}

The ST-graph for the SIR example is given in Fig.~\ref{fig:stGraph} (a). It is straightforward to check that its connected components are such that exactly one term in each $G_i$ is active in the initial combination of therapies $C = T1_{off} \ | \ T2_{off}$; and that for each $G_i$ there is no action $\pi$ whose reactants involve more than one term in $G_i$. Therefore the set of therapy definitions $T$ of the SIR example is well-formed and $\mathcal{T} = \lbrace \lbrace T1_{off}, T1_{on}\rbrace, \lbrace T2_{off}, T2_{on}\rbrace \rbrace$ is the set of switching terms.

The combination of the different switching therapies allows us to extract the possible discrete \textit{modes} in the hybrid semantics. We define a \textit{mode graph} structure as the Cartesian product of the switching terms in a ST-graph. Figure~\ref{fig:stGraph} (b) shows the mode graph for the SIR example.
\begin{defn}[Mode graph]
Let $\mathcal{M} = (S,P,T,C)$ be a D-CGF model and $\mathcal{C(G_M)} = \lbrace G_1, \ldots , G_n \rbrace$ be the set of connected components of the associated ST-graph. The mode graph $MG_\mathcal{M}$ is defined by the Cartesian product $G_1 \times \ldots \times G_n$.
\end{defn}
\begin{figure}
\centering
\subfloat[]{\includegraphics[height=3cm]{img/st-graph.pdf}}\hspace{1cm}
\subfloat[]{\includegraphics[height=3cm]{img/st_prod.pdf}}
\caption{ST-graph for the SIR example (a) and its corresponding mode graph (b). The connected components of the ST-graph are such that there is exactly one term in each $G_i$ active in $C$ and for each $G_i$, no action $\pi$ involves more than one $G_i$-term in $react(\pi)$. This ensures that $T$ is well formed and the connected nodes are switching therapies. In (b) four distinct discrete modes are generated from the Cartesian product of $G_1$ and $G_2$.}
\label{fig:stGraph}
\end{figure}
Note that transitions between modes could have associated exponentially distributed delays determined by the stochastic rates of the algebraic specification. In our context, we omit them assuming that such transitions are instantaneous, also because discrete modes will represent the control inputs of the hybrid dynamical system semantics.

As pointed out also in~\cite{bortolussi2009hybrid}, Petri Nets could have been similarly applied in order to study the conservation properties in a switching therapy. In particular, a set of switching therapies could be seen as a set of \textit{strictly conservative} (i.e.\ constant number of tokens) and \textit{1-bounded} (i.e.\ the maximum number of token is 1) Petri Nets.
%\iffalse
%The procedure for extracting the switching therapies from a D-CGF model $\mathcal{M} = (S,P,T,C)$ is given in Algorithm~\ref{algorithm:st}.
%\begin{algorithm}
%\begin{algorithmic}[1]
%\STATE Calculate the stoichiometric matrix $M$.
%\STATE Check if the conditions on $M$ at Def.~\ref{def:necessary} are met. If not, $T$ is not well-formed.
%\STATE Build the ST-graph associated to $\mathcal{M}$, $\mathcal{G_M}$
%\STATE Calculate the set of connected components of $\mathcal{G_M}$, $\mathcal{C(G_M)} = \lbrace G_1 = (T_1,E_1), \ldots , G_n = (T_n,E_n) \rbrace$.
%\STATE Check if the set $\mathcal{T}=\lbrace T_1, \ldots, T_n \rbrace$ is a partition of $T$. If not, $T$ is not well-formed.
%\STATE For each $T_i \in \mathcal{T}$, check if exactly one of the terms $U \in T_i$ is active in the initial combination of therapies: $\#(T_i, C) = 1$. If not, $T$ is not well-formed.
%\STATE For each $T_i \in \mathcal{T}$, and for each action $\pi \in \Pi$, check if the number of $T_i$-terms in $react(\pi)$ is at most one. If not, $T$ is not well-formed.
%\STATE $T$ is well-formed and each $T_i \in \mathcal{T}$ is a switching therapy.
%\end{algorithmic}
%\caption{Given a D-CGF model $\mathcal{M} = (S,P,T,C)$, calculates the set of switching therapies $\mathcal{T}$}
%\label{algorithm:st}
%\end{algorithm}
%For the sake of simplicity, we assume that condition 2 of Definition~\ref{def:st} is met, since it is quite straightforward to check if there are repeated $T$-terms in the initial combination of therapies $C$. Moreover, we do not give any implementation details on the function $computeStochMatrix(\mathcal{M},T)$ that returns the stoichiometric matrix associated to $\mathcal{M}$ and restricted to the set of terms $T$.
%\begin{algorithm}
%\begin{algorithmic}
%\STATE $\mathcal{T_M} \gets \emptyset$
%\STATE $M_{|T} \gets computeStochMatrix(\mathcal{M},T)$
%\FORALL {$\pi \in \Pi$}
%	\STATE $temp \gets \mathsf{nil}$
%	\STATE $found \gets \mathsf{false}$
%	\FORALL {$U_i \in T$}
%		\IF {$M(U_i,\pi) \neq 0$}
%			\IF {$temp = \mathsf{nil}$}
%				\STATE $temp \gets U_i$			
%			\ELSIF {$\neg found \wedge M_{|T}[U_i,\pi]+M_{|T}[temp,\pi]=0$}
%				\STATE $found \gets \mathsf{true}$
%				
%				\STATE $\mathcal{T_M} \gets \mathcal{T_M} \cup \lbrace U_i, temp\rbrace$
%			\ELSE \RETURN $\emptyset$
%		\ENDIF			
%		\ENDIF
%	\ENDFOR	
%	\IF {$temp \neq \mathsf{nil} \wedge \neg found$}
%		\RETURN $\emptyset$
%	\ENDIF
%\ENDFOR
%\end{algorithmic}
%\caption{Given a D-CGF model $\mathcal{M} = (S,P,T,C)$, calculates the set of switching therapies $\mathcal{T_M}$}
%\label{algorithm:st}
%\end{algorithm}
%\fi
\subsubsection{Hybrid dynamical system of D-CGF models}
Here we describe the semantics of a D-CGF model in terms of a particular class of hybrid dynamical systems, called \textit{controlled switched systems (CSS)}~\cite{lunze2009handbook}. While in the general formulation of hybrid dynamical system the external control input and the discrete operation mode are distinct, in a CSS the external controller produces a switching signal (i.e.\ the discrete mode) that is given in input to the plant (i.e.\ the controlled system). The basic form of a controlled switched system is the following:
\begin{align*}
\mathbf{\dot{x}} = & \mathbf{f}(\mathbf{x}, \mathbf{q})\\
\mathbf{y} = & \mathbf{g}(\mathbf{x}, \mathbf{q}),
\end{align*}
where $\mathbf{x} \in \mathbb{R}^n$ is the continuous state, $\mathbf{f}$ is a vector-valued smooth function, $\mathbf{g}$ is the output function and $\mathbf{q}$ is the discrete operation mode. The above equations are often written with the form $\mathbf{\dot{x}} = \mathbf{f_q}(\mathbf{x}) \quad \mathbf{y} = \mathbf{g_q}(\mathbf{x})$, stressing the different dynamics and outputs under different operation modes. Note that a controlled switched system can also be seen as a hybrid dynamical system with controlled discrete inputs. Figure~\ref{fig:css_control_loop} shows the control loop in a CSS.
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{img/controlled_switched.pdf}
\caption{Control loop in a controlled switched system. The plant defines the piecewise smooth dynamics of the system $\mathbf{\dot{x}}$ and the observable output $\mathbf{y}$. The output could represent patient's medical data obtained after a visit. Depending on $\mathbf{y}$, the controller acts as therapy scheduler and determines the operation mode $\mathbf{q}$ (e.g.\ a particular combination of therapies) of the plant.}
\label{fig:css_control_loop}
\end{figure}
We show how to derive a CSS from a D-CGF model, by following the method elaborated by Cardelli for the CGF algebra as regards the continuous part of the semantics, but considering also the switching therapies and the mode graph as regards  the discrete part. Recalling that each node of a mode graph $MG_\mathcal{M}=(\lbrace q_1=\lbrace U_{11}, \ldots, U_{1n}\rbrace, \ldots , q_k=\lbrace U_{k1}, \ldots, U_{km}\rbrace\rbrace, E)$ models a particular combination of therapies, a mode $q_i$ being active means that therapies $U\in q_i$ are active and that those in $T \setminus q_i$ are not. Therefore a modified version of the rate vector will be used: $$\phi_{q_i} = \phi[\faktor{\mathbf{1}}{U},\faktor{\mathbf{0}}{V}], \ U \in q_i, V \in T \setminus q_i,$$ which is obtained by zeroing all the elements of $\phi$ that contain a therapy term not belonging to $q_i$, and omitting (substitution with 1) those in $q_1$ . The resulting controlled switched system would be:
\[
\mathbf{\dot{x}} = \mathbf{f}_{q_i}(\mathbf{x}) = M_{|S} \cdot \phi_{q_i},
\]
where $q_i \in MG_\mathcal{M}$ is a node of the mode graph, and $M_{|S}$ is the stoichiometric matrix restricted to the set of species $S$. The operation modes $\mathbf{q}$, the modified rate vector $\phi_{\mathbf{q}}$ and the resulting CSS $\mathbf{\dot{x}}$ of the SIR example are shown below.
{\footnotesize
\[
\phi_{\mathbf{q}} = \bordermatrix{~ &\tau_{S_1} & \tau_{I_1} & \tau_{R_1} & \tau_{S_2} & \tau_{I_2} & \tau_{R_2} & \tau_{I_3} & \tau_{1on} & \tau_{1off} & \tau_{2on} & \tau_{2off} & i & j & h \cr
q_1& bS     & bI     & bR     & \mu S  & \mu I  & \mu R  &  \nu I  & r1_{on} & 0 & r2_{on}  &  0 & \beta SI & 0 & 0 \cr
q_2& \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &  \vdots  & 0 & r1_{off} & r2_{on} &  0 & \vdots & \rho S &  0\cr
q_3& \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &  \vdots  & r1_{on} & 0 & 0 & r2_{off} & \vdots & 0 &  k I\cr
q_4& bS     & bI     & bR     & \mu S  & \mu I  & \mu R  &  \nu I  & 0 & r1_{off} & r2_{off} & 0 & \beta SI & \rho S & k I\cr}\]\[
\mathbf{q}=\begin{bmatrix}
q_1 = (T1_{off}, T2_{off})\\
q_2 = (T1_{on}, T2_{off})\\
q_3 = (T1_{off}, T2_{on})\\
q_4 = (T1_{on}, T2_{on})
\end{bmatrix}
\qquad
\mathbf{\dot{x}}(q_1)= \left\{ 
  \begin{array}{l}
	\dot{S} = bN - \beta S I - \mu S\\
	\dot{I} = \beta S I - \mu I - \nu I\\
	\dot{R} = \nu I - \mu R\\
  \end{array} \right.
  \qquad
\mathbf{\dot{x}}(q_2)= \left\{ 
  \begin{array}{l}
	\dot{S} = bN - \beta S I - \mu S - \rho S\\
	\dot{I} = \beta S I - \mu I - \nu I\\
	\dot{R} = \nu I - \mu R + \rho S\\
  \end{array} \right.
\]\[
\mathbf{\dot{x}}(q_3)= \left\{ 
  \begin{array}{l}
	\dot{S} = bN - \beta S I - \mu S \\
	\dot{I} = \beta S I - \mu I - (\nu + k) I\\
	\dot{R} = (\nu + k) I - \mu R\\
  \end{array} \right.
    \qquad
\mathbf{\dot{x}}(q_4)= \left\{ 
  \begin{array}{l}
	\dot{S} = bN - \beta S I - \mu S - \rho S\\
	\dot{I} = \beta S I - \mu I - (\nu + k) I\\
	\dot{R} = (\nu + k) I - \mu R + \rho S\\
  \end{array} \right.
\]
}\noindent Alternatively, we could consider the equivalent formulation which exploits the sets of switching terms $\lbrace T1_{off}, T1_{on}\rbrace$ and $\lbrace T2_{off}, T2_{on}\rbrace$ as discrete control variables in the following way.
{\footnotesize
\[
T_1(q_i)= \left\{ 
  \begin{array}{ll}
	1 & \quad \text{if } T1_{on} \in q_i\\
	0 & \quad \text{if } T1_{off} \in q_i\\
  \end{array} \right.
  \qquad
  T_2(q_i)= \left\{ 
  \begin{array}{ll}
	1 & \quad \text{if } T2_{on} \in q_i\\
	0 & \quad \text{if } T2_{off} \in q_i\\
  \end{array} \right.
  \qquad
  \mathbf{\dot{x}}(q_i)= \left\{ 
  \begin{array}{l}
	\dot{S} = bN - \beta S I - \mu S - T_1(q_i)\rho S\\
	\dot{I} = \beta S I - \mu I - (\nu + T_2(q_i)k) I\\
	\dot{R} = (\nu + T_2(q_i)k) I - \mu R + T_1(q_i)\rho S\\
  \end{array}\right.
\]
}\noindent Note that this form is applicable because in the SIR example each switching term contains exactly two terms and can consequently have assigned a value in $\lbrace 0,1 \rbrace$.
\section{Hybrid semantics and optimal control of infectious processes}\label{sect:sim}
%\subsection{The SIR case study}
%In this section we report the simulation results of the SIR model with parameters that are typical of the measles~\cite{stone2000theoretical} and define the optimal control law for therapy scheduling by means of Model Predictive Control (MPC)~\cite{bemporad1999control}. A similar problem has been addressed in~\cite{suzuki2010piecewise} for the scheduling of hormone therapy in a model for prostate cancer.
In this section we define the optimal control law for therapy scheduling in the SIR model by means of Model Predictive Control (MPC)~\cite{bemporad1999control}, and we report the optimally controlled solutions with parameters that are typical of the measles~\cite{stone2000theoretical}. A similar problem has been addressed in~\cite{suzuki2010piecewise} for the scheduling of hormone therapy in a model for prostate cancer.
%
%
%Figure~\ref{fig:SIR_simulation} shows the numerical results of the ODE model obtained by running the four different therapy administration strategies (the discrete modes), introduced in Sect.~\ref{subsect:hybridsemantics}. Under therapy $T_2$, which increases the recovery rate, larger oscillations are observed w.r.t. the scenario without therapies and the dynamics unexpectedly lead to a stable state with a higher percentage of susceptible population and a lower percentage of recovered. On the other hand, by administering therapy $T_1$ stability is reached before and higher percentages of recovered individuals are recorded. Finally, the combination of $T_1$ and $T_2$ gives the same outcome as in the scenario with only therapy $T_1$ active.
%\begin{figure}
%\centering
%\includegraphics[width=\textwidth]{img/sir_plots.pdf}
%\caption{Simulation results of the SIR model under four different therapy administration strategies. Parameters are: $S(0)=0.9, \ I(0)=0.1, \ R(0)=0,\ b = \mu = 0.02, \ \beta=1800, \ \nu = 100, \ \rho=0.5, \ k=50$.}
%\label{fig:SIR_simulation}
%\end{figure}

In order to cope with the nonlinear dynamics of the SIR model, the optimal scheduling of multiple therapies is computed by solving an on-line MPC problem in the Multi-Parametric Toolbox (MPT)~\cite{mpt}. The optimal control law can be formulated by means of the following constrained finite-time optimal control (CFTOC) problem, where the original continuous-time dynamics is simply discretized with the Euler method:
{\footnotesize
\[\begin{array}{ll}
%&\underset{\mathbf{u}}{min} \Vert P\mathbf{x}(T) \Vert_{1} + \sum_{k=0}^{T-1} \Vert R\mathbf{u}(k) \Vert_{1} + \Vert Q\mathbf{x}(k) \Vert_{1}\\
&\underset{\mathbf{u}}{min} \sum_{k=0}^{T-1} \Vert R\mathbf{u}(k) \Vert_{1} + \Vert Q\mathbf{x}(k) \Vert_{1}\\
subj. \ to \quad & \mathbf{x}(k) \in [0,1]^3\\
& \mathbf{u}(k) \in \lbrace 0,1\rbrace^2\\
& \mathbf{x}(T) \in \mathcal{X}_{set}\\
&\mathbf{x}(k+\Delta t)= A\mathbf{x}(k) + B\mathbf{u}(k)\\
&\forall k = 0,\ldots , T
\end{array}
\quad
\begin{array}{ll}
A = &\begin{bmatrix}
1 + (b -\beta I(k)-\mu)\Delta t & b\Delta t & b\Delta t\\
0 & 1 + (\beta S(k)-\mu - \nu)\Delta t & 0\\
0 & \nu\Delta t & 1 -\mu\Delta t
\end{bmatrix}\\
B = &\begin{bmatrix}
-\rho S(k) \Delta t & 0\\
0 & -k I(k) \Delta t\\
\rho S(k) \Delta t & k I(k) \Delta t\end{bmatrix} \quad 
Q = \begin{bmatrix} 1 & 0 & 0\\ 0 & 10 & 0\\ 0 &0&0.5\end{bmatrix}
\end{array}
\]
}\noindent where $\Delta t = 1/365$ is the discrete time step corresponding to one day; $\mathbf{u}(k)= [T1(k) \ T2(k)]^T \in \lbrace 0,1\rbrace^2$ is the discrete input (the combination of the two therapies) and $\mathbf{x}(k)= [S(k) \ I(k) \ R(k)]^T \in [0,1]^3$ is the systems' continuous state vector at time $k$; $T=3$ is the prediction horizon (corresponding to three days); $\mathcal{X}_{set}$ is a polytope identifying the set of terminal states; $\Vert \mathbf{x} \Vert_{1} = |S| + |I| + |R|$ and $\Vert \mathbf{u} \Vert_{1} = |T1| + |T2|$; $R$ and $Q$ are the weights on the controlled inputs and on the system state, resp. The terminal set $\mathcal{X}_{set}$ represents an invariant on the final states and can be described by the convex hull of its vertices: $\mathcal{X}_{set} = Conv(\lbrace [1 \ 0 \ 0]^T, [0 \ 0 \ 1]^T\rbrace)$. In this way, the terminal set contains all the states $[S \ I \ R]^T$ such that $I=0$ (no infected individuals), $S,R \geq 0$ and $S+R = 1$ (we do not exceed the population limit). Matrix $Q$ has been set so that the state variable Infected has associated a weight equals to $10$ and thus is slightly penalized over variables Susceptible (weight $1$) and Recovered (weight $0.5$). Finally, we assume for simplicity that the internal state can be observed and we set the output vector $\mathbf{y}$ to the state vector $\mathbf{x}$.

Matrix $R$ allows us to implement different therapy scheduling strategies. The higher $R$, the higher penalty is given to drug dosage. In particular if $R$ is assigned a large value, then the strategy is to avoid drug dosage as much as possible. This can be the case of a therapy with severe side effects or of a patient in an early stage of the disease that can be treated even with a small dosage. On the contrary, a low value to $R$ leads to a strategy where drug dosage is much more prominent (e.g.\ little side effects or mature stage of disease). Additionally we can assign different weights to different drugs according to their therapeutic impact on the patient, thus possibly enabling a wider spectrum of control strategies. We computed the optimally controlled solutions considering three scenarios obtained by varying the weights on the controlled therapies:
\begin{enumerate}
\item Low penalties to $T1$ and $T2$: $R=\begin{bmatrix} 0.1 & 0 \\ 0 & 0.1\end{bmatrix}$
\item High penalty to $T1$: $R=\begin{bmatrix} 100 & 0 \\ 0 & 0.1\end{bmatrix}$
\item High penalty to $T2$: $R=\begin{bmatrix} 0.1 & 0 \\ 0 & 100\end{bmatrix}$
\end{enumerate}

Results reported in Figure~\ref{fig:SIR_MPT} tell that solutions under scenario 1 and scenario 2 are the same (Fig.~\ref{fig:SIR_MPT} (a)), suggesting that therapy $T1$ does not need to be administered in the optimal strategy of treatment in order to fulfil the terminal constraints. Indeed Fig.~\ref{fig:SIR_MPT} (b) shows that when $T2$ has a high weight (scenario 3), no control moves are performed, so indicating that $T1$ is ineffective with respect to $T2$.
\begin{figure}
\centering
\subfloat[Scenarios 1,2]{\includegraphics[width=0.3\textwidth, height=0.5\textwidth]{img/s1_global.pdf}}\hspace*{1cm}
\subfloat[Scenario 3]{\includegraphics[width=0.3\textwidth, height=0.5\textwidth]{img/s3_local.pdf}}
\caption{Optimally controlled solutions for the SIR model with therapies as controlled discrete inputs simulated for 15 days. Parameters are: $S(0)=0.3, \ I(0)=0.7, \ R(0)=0,\ b = \mu = 0.02, \ \beta=1800, \ \nu = 100, \ \rho=0.5, \ k=50$. Results have been computed by using a global optimal nonlinear solver.}
\label{fig:SIR_MPT}
\end{figure}
%\begin{figure}
%\centering
%\begin{tabular}{m{0.2\textwidth}|m{0.35\textwidth}m{0.35\textwidth}}
%&Local solver& Global solver\\
%\hline
%$\begin{array}{cc}
%R=\begin{bmatrix} 0.1 & 0 \\ 0 & 0.1\end{bmatrix} \\
%\\
%R=\begin{bmatrix} 100 & 0 \\ 0 & 0.1\end{bmatrix}
%\end{array}$&\subfloat[]{\includegraphics[width=0.3\textwidth, height=0.5\textwidth]{img/s1_local.pdf}}&
%\subfloat[]{\includegraphics[width=0.3\textwidth, height=0.5\textwidth]{img/s1_global.pdf}}\\\hline
%$R=\begin{bmatrix} 0.1 & 0 \\ 0 & 100\end{bmatrix}$ &\subfloat[]{\includegraphics[width=0.3\textwidth, height=0.5\textwidth]{img/s3_local.pdf}}&
%\subfloat[]{\includegraphics[width=0.3\textwidth, height=0.5\textwidth]{img/s3_local.pdf}}
%\end{tabular}
%\caption{}
%\label{fig:SIR_MPT}
%\end{figure}

\subsection{Osteomyelitis, a truly multiscale infection}
Osteomyelitis is a a bone pathology caused by bacteria infection (mostly \textit{Staphylococcus aureus}) that alters the bone remodelling process and that rapidly leads to severe bone loss, necrosis of the affected portion, and it may even spread to other parts of the body. Therefore it is one of the most appropriate examples of multiscale infection, since it affects multiple biological scales. Indeed infection by \textit{S. aureus} starts at the intracellular level, but it subsequently involves also the tissue and the organ level. \textit{Bone remodelling (BR)} is the cellular-level process by which the bone is continuously renewed as the result of an alternation of bone resorption, conducted by cells called {\it osteoclasts}, and bone formation, conducted by {\it osteoblasts}. They form together the so-called \textit{Basic Multicellular Unit (BMU)}, i.e.\ an ensemble of osteoclasts and osteoblasts that dissolve an area of the bone surface and then fill it. The BMU could be considered as an emerging scale intermediate between the cell and tissue levels. Osteomyelitis induces a severe inflammatory response followed by progressive bone destruction and loss of the vasculature and with a persistent chronic infection; this is further complicated by the rapid emergence of resistant strains of \textit{S. aureus}. It has been shown that the infection prevents proliferation, induces apoptosis and inhibits mineralisation of cultured osteoblasts. Although effective treatment of this disease is very difficult, one of most used drug is the fusidic acid that acts as a bacteriostatic agent, and is usually combined with other antibiotics. Based on a recent work of the authors about the modelling of osteomyelitis and the comparison of different treatments~\cite{lio2012osteomyelitis}, we can define a hybrid version of that model and a treatment strategy combining antibiotic and anti-inflammatory therapies. The model describes how the bone remodelling dynamics of osteoblasts' ($Ob$) and osteoclasts' ($Oc$) population is affected by the \textit{S.aureus} ($B$).
{\small
\[
\mathbf{x}=\begin{bmatrix}
O_{c}\\O_{b}\\B
\end{bmatrix}\quad
\mathbf{u}=\begin{bmatrix}
T1\\T2
\end{bmatrix} \quad
\begin{matrix}
\dot{O_{c}}=\alpha_{1}O_{c}^{g_{11}(1+f_{11}\frac{B}{s})}O_{b}^{g_{21}(1 +T2\cdot k_i- f_{21}\frac{B}{s})}-\beta_{1}O_{c},\\
\dot{O_{b}}=\alpha_{2}O_{c}^{g_{12}/(1+f_{12}\frac{B}{s})}O_{b}^{g_{22}-f_{22}\frac{B}{s}}- \beta_{2} O_{b},\\
\dot{B}=(1-T1)\gamma_{B}B\cdot log(\frac{s}{B})
\quad
y = -k_1\cdot Oc + k_2 \cdot Ob,
\end{matrix}
\]
}\noindent where $\alpha_{i}$ and $\beta_{i}$ are growth and death rates; $g_{ij}$ describes the effectiveness of autocrine and paracrine regulation, i.e.\ the chemical interactions between osteoblasts and osteoclasts; $f_{ij}$ models the impact of the infection on the autocrine and paracrine regulation; $T1$ is the discrete input modelling the antibiotic dosage; $T2$ models the dosage of a anti-inflammatory therapy; $k_i$ is the 
effectiveness of therapy $T2$; the bacterial population follows a Gompertz curve with carrying capacity $s$ and growth rate $\gamma_{B}$; and the output function $y$ represents the bone density calculated as a function of $Oc$, $Ob$ and their resorption ($k_1$) and formation ($k_2$) rates. Bone Mineral Density (BMD) measurement are often taken in medical practice, thus $y$ is the quantity observed by the controller (a doctor who administers the treatment) after each visit.

The optimal scheduling of bacteriostatic and anti-inflammatory therapies could be formulated with a control law similar to the SIR example. Since the combination of antibiotic and anti-inflammatory drugs has never been applied in medical practice, we could set the weight matrix of input therapies $R$ so that $R_{12}$ and $R_{21}$ are much larger that $R_{11}$ and $R_{22}$. In this way, we give a high penalty to the combination of therapies $T1$ and $T2$ w.r.t. $T1$ and $T2$ taken individually.

\section{Discussion and conclusion}\label{sect:conclusions}
The field of predictive models in biomedicine is challenged by the need of a better comprehension of the phenomenon of scales in the biological organisation, particularly their role in the transition between health and disease conditions. The multiscale modelling of molecules-cell-tissue-organ-body interactions is a key step in the process of identifying the most important parameters acting in a disease state and their calibration, linking basic research and clinical practice (therapies). Here we discuss how a description based on hybrid dynamical systems is beneficial to the formulation of a multiscale modelling in biomedical processes. A second instance of the utility of hybrid dynamics approach stems from the presence of multiple controls in biological systems. These controls could be framed as occurring naturally (the immune system response) or induced after the administration of a therapy. It is very often the case of multiple treatments, i.e. switching between therapies or combination of therapies which is the case considered in this study. In particular, we show the suitability of the D-CGF high-level process algebra to accomplish the task: here, its semantics is given in terms of a hybrid dynamical system where the large number of cell populations provides the basis for a continuous modelling approach and the dosage of multiple therapies is implemented as discrete switches.
%At the lowest scale, the control is implemented on top of a Gompertz curve which is known from experiments reported in literature to provide a good description of the developing of bacterial infections (here staphylococcus aureus, very common in hospitals).

We believe that this approach opens several interesting directions in basic and clinical research.
%which will constitute our work for the next months.
First, it is general and could be used for different tissues and organs or multiorgan diseases. Second, it is possible to extend the mathematical formulation to all the scales of biological organisation involved in the infection and recovery conditions. Third, the therapy could be complex and built on several nested and hierarchical protocols.
One final comment: given the richness of examples provided by biological processes and medical therapies, there is basis for deriving interesting theories. We remind the \textit{bon mot} of Stan Ulam: Ask not what mathematics can do for biology. Ask what biology can do for Mathematics.

%\fi
\bibliographystyle{eptcs}
\bibliography{hsb}
\end{document}

hsb_camera_ready_ver2.bbl
\begin{thebibliography}{10}
\providecommand{\bibitemdeclare}[2]{}
\providecommand{\surnamestart}{}
\providecommand{\surnameend}{}
\providecommand{\urlprefix}{Available at }
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\href}[2]{\texttt{#2}}
\providecommand{\urlalt}[2]{\href{#1}{#2}}
\providecommand{\doi}[1]{doi:\urlalt{http://dx.doi.org/#1}{#1}}
\providecommand{\bibinfo}[2]{#2}

\bibitemdeclare{article}{bagnoli2007risk}
\bibitem{bagnoli2007risk}
\bibinfo{author}{F.~\surnamestart Bagnoli\surnameend},
  \bibinfo{author}{P.~\surnamestart Li\`o\surnameend} \&
  \bibinfo{author}{L.~\surnamestart Sguanci\surnameend} (\bibinfo{year}{2007}):
  \emph{\bibinfo{title}{Risk perception in epidemic modeling}}.
\newblock {\sl \bibinfo{journal}{Physical Review E}}
  \bibinfo{volume}{76}(\bibinfo{number}{6}), p. \bibinfo{pages}{061904},
  \doi{10.1103/PhysRevE.76.061904}.

\bibitemdeclare{article}{bartocci2012tcsb}
\bibitem{bartocci2012tcsb}
\bibinfo{author}{E.~\surnamestart Bartocci\surnameend},
  \bibinfo{author}{P.~\surnamestart Li\`o\surnameend},
  \bibinfo{author}{M.~\surnamestart Merelli\surnameend} \&
  \bibinfo{author}{N.~\surnamestart Paoletti\surnameend}
  (\bibinfo{year}{2012}): \emph{\bibinfo{title}{Multiple Verification in
  Complex Biological Systems: the Bone Remodelling Case Study}}.
\newblock {\sl \bibinfo{journal}{Transactions on Computational Systems Biology.
  To appear}}.

\bibitemdeclare{article}{batt2005qualitative}
\bibitem{batt2005qualitative}
\bibinfo{author}{G.~\surnamestart Batt\surnameend},
  \bibinfo{author}{D.~\surnamestart Ropers\surnameend},
  \bibinfo{author}{H.~\surnamestart De~Jong\surnameend},
  \bibinfo{author}{J.~\surnamestart Geiselmann\surnameend},
  \bibinfo{author}{M.~\surnamestart Page\surnameend} \&
  \bibinfo{author}{D.~\surnamestart Schneider\surnameend}
  (\bibinfo{year}{2005}): \emph{\bibinfo{title}{Qualitative analysis and
  verification of hybrid models of genetic regulatory networks: Nutritional
  stress response in Escherichia coli}}.
\newblock {\sl \bibinfo{journal}{Hybrid Systems: Computation and Control}}, pp.
  \bibinfo{pages}{134--150}, \doi{10.1007/978-3-540-31954-2\_9}.

\bibitemdeclare{article}{bemporad1999control}
\bibitem{bemporad1999control}
\bibinfo{author}{A.~\surnamestart Bemporad\surnameend} \&
  \bibinfo{author}{M.~\surnamestart Morari\surnameend} (\bibinfo{year}{1999}):
  \emph{\bibinfo{title}{Control of systems integrating logic, dynamics, and
  constraints}}.
\newblock {\sl \bibinfo{journal}{Automatica}} \bibinfo{volume}{35}, pp.
  \bibinfo{pages}{407--428}, \doi{10.1016/S0005-1098(98)00178-2}.

\bibitemdeclare{inproceedings}{bortolussi2010hybrid}
\bibitem{bortolussi2010hybrid}
\bibinfo{author}{L.~\surnamestart Bortolussi\surnameend},
  \bibinfo{author}{V.~\surnamestart Galpin\surnameend},
  \bibinfo{author}{J.~\surnamestart Hillston\surnameend} \&
  \bibinfo{author}{M.~\surnamestart Tribastone\surnameend}
  (\bibinfo{year}{2010}): \emph{\bibinfo{title}{Hybrid semantics for PEPA}}.
\newblock In: {\sl \bibinfo{booktitle}{Quantitative Evaluation of Systems
  (QEST), 2010 Seventh International Conference on the}},
  \bibinfo{organization}{IEEE}, pp. \bibinfo{pages}{181--190},
  \doi{10.1109/QEST.2010.31}.

\bibitemdeclare{article}{bortolussi2009hybrid}
\bibitem{bortolussi2009hybrid}
\bibinfo{author}{L.~\surnamestart Bortolussi\surnameend} \&
  \bibinfo{author}{A.~\surnamestart Policriti\surnameend}
  (\bibinfo{year}{2009}): \emph{\bibinfo{title}{Hybrid Dynamics of Stochastic
  $\pi$-calculus}}.
\newblock {\sl \bibinfo{journal}{Mathematics in Computer Science}}
  \bibinfo{volume}{2}(\bibinfo{number}{3}), pp. \bibinfo{pages}{465--491},
  \doi{10.1007/s11786-008-0065-3}.

\bibitemdeclare{inproceedings}{cardelli2008processes}
\bibitem{cardelli2008processes}
\bibinfo{author}{L.~\surnamestart Cardelli\surnameend} (\bibinfo{year}{2008}):
  \emph{\bibinfo{title}{From processes to odes by chemistry}}.
\newblock In: {\sl \bibinfo{booktitle}{Fifth Ifip International Conference On
  Theoretical Computer Science--Tcs 2008}}, \bibinfo{organization}{Springer},
  pp. \bibinfo{pages}{261--281}, \doi{10.1007/978-0-387-09680-3\_18}.

\bibitemdeclare{article}{cardelli2008computational}
\bibitem{cardelli2008computational}
\bibinfo{author}{L.~\surnamestart Cardelli\surnameend} \&
  \bibinfo{author}{G.~\surnamestart Zavattaro\surnameend}
  (\bibinfo{year}{2008}): \emph{\bibinfo{title}{On the computational power of
  biochemistry}}.
\newblock {\sl \bibinfo{journal}{Algebraic Biology}}, pp.
  \bibinfo{pages}{65--80}, \doi{10.1007/978-3-540-85101-1\_6}.

\bibitemdeclare{article}{galpin2009hype}
\bibitem{galpin2009hype}
\bibinfo{author}{V.~\surnamestart Galpin\surnameend},
  \bibinfo{author}{L.~\surnamestart Bortolussi\surnameend} \&
  \bibinfo{author}{J.~\surnamestart Hillston\surnameend}
  (\bibinfo{year}{2009}): \emph{\bibinfo{title}{HYPE: a process algebra for
  compositional flows and emergent behaviour}}.
\newblock {\sl \bibinfo{journal}{CONCUR 2009-Concurrency Theory}}, pp.
  \bibinfo{pages}{305--320}, \doi{10.1007/978-3-642-04081-8\_21}.

\bibitemdeclare{inproceedings}{grosu2011cardiac}
\bibitem{grosu2011cardiac}
\bibinfo{author}{R.~\surnamestart Grosu\surnameend},
  \bibinfo{author}{G.~\surnamestart Batt\surnameend},
  \bibinfo{author}{F.~\surnamestart Fenton\surnameend},
  \bibinfo{author}{J.~\surnamestart Glimm\surnameend},
  \bibinfo{author}{C.~\surnamestart Le~Guernic\surnameend},
  \bibinfo{author}{S.~\surnamestart Smolka\surnameend} \&
  \bibinfo{author}{E.~\surnamestart Bartocci\surnameend}
  (\bibinfo{year}{2011}): \emph{\bibinfo{title}{From cardiac cells to genetic
  regulatory networks}}.
\newblock In: {\sl \bibinfo{booktitle}{Computer Aided Verification}},
  \bibinfo{organization}{Springer}, pp. \bibinfo{pages}{396--411},
  \doi{10.1007/978-3-642-22110-1\_31}.

\bibitemdeclare{inproceedings}{henzinger1996theory}
\bibitem{henzinger1996theory}
\bibinfo{author}{T.A. \surnamestart Henzinger\surnameend}
  (\bibinfo{year}{1996}): \emph{\bibinfo{title}{The theory of hybrid
  automata}}.
\newblock In: {\sl \bibinfo{booktitle}{Logic in Computer Science, 1996.
  LICS'96. Proceedings., Eleventh Annual IEEE Symposium on}},
  \bibinfo{organization}{IEEE}, pp. \bibinfo{pages}{278--292},
  \doi{10.1109/LICS.1996.561342}.

\bibitemdeclare{article}{kermack1932contributions}
\bibitem{kermack1932contributions}
\bibinfo{author}{W.O. \surnamestart Kermack\surnameend} \&
  \bibinfo{author}{A.G. \surnamestart McKendrick\surnameend}
  (\bibinfo{year}{1932}): \emph{\bibinfo{title}{Contributions to the
  mathematical theory of epidemics. II. The problem of endemicity}}.
\newblock {\sl \bibinfo{journal}{Proceedings of the Royal society of London.
  Series A}} \bibinfo{volume}{138}(\bibinfo{number}{834}), pp.
  \bibinfo{pages}{55--83}, \doi{10.1098/rspa.1932.0171}.

\bibitemdeclare{misc}{mpt}
\bibitem{mpt}
\bibinfo{author}{M.~\surnamestart Kvasnica\surnameend},
  \bibinfo{author}{P.~\surnamestart Grieder\surnameend} \&
  \bibinfo{author}{M.~\surnamestart Baoti\'{c}\surnameend}
  (\bibinfo{year}{2004}): \emph{\bibinfo{title}{{Multi-Parametric Toolbox
  (MPT)}}}.
\newblock \urlprefix\url{http://control.ee.ethz.ch/~mpt/}.

\bibitemdeclare{inproceedings}{lapin2011shave}
\bibitem{lapin2011shave}
\bibinfo{author}{M.~\surnamestart Lapin\surnameend},
  \bibinfo{author}{L.~\surnamestart Mikeev\surnameend} \&
  \bibinfo{author}{V.~\surnamestart Wolf\surnameend} (\bibinfo{year}{2011}):
  \emph{\bibinfo{title}{SHAVE: stochastic hybrid analysis of markov population
  models}}.
\newblock In: {\sl \bibinfo{booktitle}{Proceedings of the 14th international
  conference on Hybrid systems: computation and control}},
  \bibinfo{organization}{ACM}, pp. \bibinfo{pages}{311--312},
  \doi{10.1145/1967701.1967746}.

\bibitemdeclare{article}{lio2012osteomyelitis}
\bibitem{lio2012osteomyelitis}
\bibinfo{author}{P.~\surnamestart Li\`o\surnameend},
  \bibinfo{author}{N.~\surnamestart Paoletti\surnameend}, \bibinfo{author}{M.A.
  \surnamestart Moni\surnameend}, \bibinfo{author}{K.~\surnamestart
  Atwell\surnameend}, \bibinfo{author}{M.~\surnamestart Merelli\surnameend} \&
  \bibinfo{author}{Viceconti \surnamestart M.\surnameend}
  (\bibinfo{year}{2012}): \emph{\bibinfo{title}{{Modelling osteomyelitis}}}.
\newblock {\sl \bibinfo{journal}{BMC Bioinformatics, supplement on NETTAB 2011
  workshop. To appear}}.

\bibitemdeclare{book}{lunze2009handbook}
\bibitem{lunze2009handbook}
\bibinfo{author}{J.~\surnamestart Lunze\surnameend},
  \bibinfo{author}{F.~\surnamestart Lamnabhi-Lagarrigue\surnameend} et~al.
  (\bibinfo{year}{2009}): \emph{\bibinfo{title}{Handbook of hybrid systems
  control: theory, tools, applications}}.
\newblock \bibinfo{publisher}{Cambridge University Press}.

\bibitemdeclare{article}{norman2003developing}
\bibitem{norman2003developing}
\bibinfo{author}{R.~\surnamestart Norman\surnameend} \&
  \bibinfo{author}{C.~\surnamestart Shankland\surnameend}
  (\bibinfo{year}{2003}): \emph{\bibinfo{title}{Developing the use of process
  algebra in the derivation and analysis of mathematical models of infectious
  disease}}.
\newblock {\sl \bibinfo{journal}{Computer Aided Systems Theory-EUROCAST 2003}},
  pp. \bibinfo{pages}{404--414}, \doi{10.1007/978-3-540-45210-2\_37}.

\bibitemdeclare{inproceedings}{paoletti2011}
\bibitem{paoletti2011}
\bibinfo{author}{N.~\surnamestart Paoletti\surnameend},
  \bibinfo{author}{P.~\surnamestart Li\`{o}\surnameend},
  \bibinfo{author}{E.~\surnamestart Merelli\surnameend} \&
  \bibinfo{author}{M.~\surnamestart Viceconti\surnameend}
  (\bibinfo{year}{2011}): \emph{\bibinfo{title}{Osteoporosis: a multiscale
  modeling viewpoint}}.
\newblock In: {\sl \bibinfo{booktitle}{Proceedings of the 9th International
  Conference on Computational Methods in Systems Biology (CMSB '11)}}, pp.
  \bibinfo{pages}{183--193}, \doi{10.1145/2037509.2037536}.

\bibitemdeclare{article}{paoletti2011spec}
\bibitem{paoletti2011spec}
\bibinfo{author}{N.~\surnamestart Paoletti\surnameend},
  \bibinfo{author}{P.~\surnamestart Li\`{o}\surnameend},
  \bibinfo{author}{E.~\surnamestart Merelli\surnameend} \&
  \bibinfo{author}{M.~\surnamestart Viceconti\surnameend}
  (\bibinfo{year}{2012}): \emph{\bibinfo{title}{Multi-level Computational
  Modeling and Quantitative Analysis of Bone Remodeling}}.
\newblock {\sl \bibinfo{journal}{IEEE/ACM Transactions on Computational Biology
  and Bioinformatics}} \bibinfo{volume}{99}(\bibinfo{number}{PrePrints}),
  \doi{10.1109/TCBB.2012.51}.

\bibitemdeclare{article}{stone2000theoretical}
\bibitem{stone2000theoretical}
\bibinfo{author}{L.~\surnamestart Stone\surnameend},
  \bibinfo{author}{B.~\surnamestart Shulgin\surnameend} \&
  \bibinfo{author}{Z.~\surnamestart Agur\surnameend} (\bibinfo{year}{2000}):
  \emph{\bibinfo{title}{Theoretical examination of the pulse vaccination policy
  in the SIR epidemic model}}.
\newblock {\sl \bibinfo{journal}{Mathematical and Computer Modelling}}
  \bibinfo{volume}{31}(\bibinfo{number}{4}), pp. \bibinfo{pages}{207--215},
  \doi{10.1016/S0895-7177(00)00040-6}.

\bibitemdeclare{article}{suzuki2010piecewise}
\bibitem{suzuki2010piecewise}
\bibinfo{author}{T.~\surnamestart Suzuki\surnameend},
  \bibinfo{author}{N.~\surnamestart Bruchovsky\surnameend} \&
  \bibinfo{author}{K.~\surnamestart Aihara\surnameend} (\bibinfo{year}{2010}):
  \emph{\bibinfo{title}{Piecewise affine systems modelling for optimizing
  hormone therapy of prostate cancer}}.
\newblock {\sl \bibinfo{journal}{Philosophical Transactions of the Royal
  Society A: Mathematical, Physical and Engineering Sciences}}
  \bibinfo{volume}{368}(\bibinfo{number}{1930}), pp.
  \bibinfo{pages}{5045--5059}, \doi{10.1098/rsta.2010.0220}.

\bibitemdeclare{book}{viceconti2011multiscale}
\bibitem{viceconti2011multiscale}
\bibinfo{author}{M.~\surnamestart Viceconti\surnameend} (\bibinfo{year}{2011}):
  \emph{\bibinfo{title}{Multiscale Modeling of the Skeletal System}}.
\newblock \bibinfo{publisher}{Cambridge University Press},
  \doi{10.1017/CBO9781139049627}.

\end{thebibliography}
