Modelling and performance analysis of clinical pathways using the stochastic process algebra PEPA

Motivation

The stroke clinical pathway investigated in this paper is for illustrative purposes. Stroke is a complex disease requiring a systematic integration of services, e.g. primary care, ambulance services, acute treatment and rehabilitation, post-acute rehabilitation, and often long-term health and care support in the community. Due to its complexity and requirement of various services, the stroke clinical pathway presents a significant challenge to existing generally fragmented health and social care services. Moving from the current fragmented approach to an integrated stroke system is a complex task. As simulation models are useful during the planning stage of complex service re-configuration, modelling methods are in demand to simulate and optimize the stroke clinical pathway. Figure 1 displays the abstraction of a general stroke clinical pathway, obtained from Charing Cross hospital of Imperial College London. A new 999 call from a patient or other hospitals initiates an instance of this pathway. The information flow is summarized as follows:

Related work

The general impact of information technology on the clinical pathway like the influence of IT support on the satisfaction of both patients and medical staff is studied [8–10]. The clinical pathway can be implemented manually on paper [11] or electronically [12]. With a good understanding of the clinical pathway, a set of evidence-based recommendations forming guidelines for clinical practice can be developed to both formalize and optimize the care process [13]. Thus IT techniques embedded into clinical pathways can efficiently decrease undesired practice variability and improve clinician performance.

Some researchers further propose methods of modelling the clinical pathway. By using quantitative models, the influence of individual resources on overall pathway performance can be evaluated directly. However, little work has been done in the area of modelling clinical pathways with the aim of quantitatively improving system performance. Although in [19], a pathway model to facilitate the stroke care planning process is introduced, without a mathematical representation, the model performance cannot be shown explicitly. The work has been done in [20] applies the workflow graph to model the clinical pathway. Using this model, a similarity function is proposed to evaluate the temporal equivalence of two clinical pathways in order to reduce a complex pathway scenario into a simple one. In [21], ontology is used to describe the clinical pathway. This model allows the hierarchical representation of the clinical pathway, e.g., a clinical pathway can be described as a combination of a high-level outcome flow and a detailed workflow with care time constraints. However, the work done in [20] and [21] cannot build a model to explicitly measure the performance of the clinical pathway, such as throughput and passage time. Some researchers focus on analysing time and resource information in clinical pathway using stochastic Petri net. In [22], Rui et al. introduce the Probabilistic Time Constraint WorkFlow Nets (PTCWF-nets) to model a process. A static analysis method is then proposed to analyse each activity's probability of meeting its time constraints (e.g., deadline) in a process before the process is actually executed. In addition, they develop a dynamic method during the execution of the process [23]. This method can update remaining activities' probabilities of successful execution whenever some activities are completed and their actual durations are known. Time schedulability is only one aspect of performance analysis in clinical pathway and it is closely related to another aspect: resource analysis. In [24], they further apply PTCWF-nets to manage resources. Their approach can schedule a clinical pathway among multiple available resources and allocate each activity to an optimal resource, i.e., the resource that has the highest probability to finish the activity in time. Authors in [25] also attempt to model clinical pathway using Petri net. They introduce performance trees to provide a standard unifying framework for expressing performance measures and performance requirements. Benefit from performance trees, the Petri net model for clinical pathway can provide estimation of steady state distribution and passage time.

Besides the work has been done to model the clinical pathway using Petri net, process algebra is also introduced to formally specify the interactions of different entities in the clinical pathway [26]. Process algebra is a mathematical framework used to describe a complex parallel system. In this framework, both the behaviour and properties of the system are described in the form of algebra, facilitating accurate definition and rigorous reasoning about the system in mathematics. PEPA is an enhanced process algebra mainly used to describe and analyse the performance of concurrent systems [14]. It inherits most characteristics of process calculus while incorporating features to specify a stochastic model, which potentially behaves as a continuous time Markov process. Comparing PEPA with other frequently used modelling tools, a queuing network offers compositionality but lacks formal definition, while a Petri net has formal definition without good compositionality. Thus in paper [27], PEPA is used to model the healthcare system. Based on PEPA, the execution duration and throughput of the clinical pathway can be evaluated. However, the role of resources, which constrain the activities and further limit pathway outcome, is not explicitly shown in this model. Therefore, this paper proposes a general clinical pathway model based on the CPP to analyse the performance of system and optimize pathway output.

The remainder of this paper is organised as follows: the Method Section introduces the architecture of ICPA, defines the CPP model, and presents key theories in conducting performance analysis on clinical pathway; the Results Section reports our experimental tests of the CPP model's effectiveness; finally, the Conclusion Section summarises our work.

Definition of CPP

Inheriting and developed from PEPA, the basic elements of CPP are state components, resource components and activities. A state component represents the status of a patient by showing in which department the patient is being treated. A resource component is used to specify the state of each resource, which can be either busy or idle. The third component is the activity, which decides transitions between different state components. As each activity relies on one or more resources, the state transition can only take place when the associated resource components are in their idle states. Therefore, the number of resources constrains the frequency of state transitions and limits the system throughput. The model defined by CPP which represents states and resources separately is suitable for resource optimization and status monitoring of patients. Definitions of CPP are as follows.

Definition 1. In CPP, the clinical pathway is represented by a five-tuple, < S, R, Act, C, F _C >, where:

Definition 2. In CPP, the set S and R are basic components. Given that these components can be commonly denoted as P or Q, the syntax of the terms in CPP can be defined as follows:

(1)

For example, $P\stackrel{def}{=}\left(\alpha ,rate\right).Q$ means that the component P becomes Q with the completion of the activity (α, rate). The expression P + Q represents that the system can behave either as P or Q. It enables all the activities of P and Q, and the first completed activity determines how the system behaves. The cooperation operator in the expression forms the basis of composition and can specify two components working cooperatively with shared activities defined in L.

Three key parts of the CPP model

Key techniques of performance analysis using CPP

Rooted in a continuous time Markov process [14], CPP model can be used to estimate the performance of clinical pathway including throughput (the number of patients that the healthcare system can serve every day) and resource utilisation (the percentage of time that a resource is in use).

The underlying stochastic process

The basic idea of a Markov process is that the distribution of time until the next state change is independent of the time which has elapsed since the last state change. Let the state transition rate between states S _i and S _j be

$q_{S_i,Sj}=\sum _{{}_{\alpha }\in S_i\underrightarrow{\alpha }{S}_j}rat{e}_{\alpha }$

(10)

If there are no direct connections between these two states, then the transition rate $q_{S_i,S_j}=0$. Let X(t) = S _i indicate that at time instance t, the system behaves as S _i . After a tiny amount of time δt, the probability that the system is in state S _j is

$Pr\left(X\left(t+\delta t\right)=Sj|X\left(t\right)=S_i\right)=q_{S_i,Sj}*\delta t+O\left(\delta t\right),i\ne j$

(11)

where O(δt) goes to zero faster than δt. Suppose that the infinitesimal generator for the process Q is a square matrix whose off-diagonal elements are $q_{S_i,S_j}$ and diagonal elements are formed as the negative sum of the non-diagonal elements of each row, $q_{S_i,S_i}=-\sum _{j\ne i}{q}_{S_i,S_j}$. The evolution of the continuous time Markov process is represented by a first order differential equation [29]:

$P^{\prime }\left(t\right)=Q*P\left(t\right)=P\left(t\right)*Q,where\mathit{ P}\left(0\right)=I$

(12)

where P(t) is a square matrix with (i, j)th entry p_{i, j}standing for the probability Pr(X(t) = S _j |X(0) = S _i ). As performance analysis is usually concerned with system behaviour over a significant period of time, it is necessary to study the steady state behaviour of the system. If the steady state distribution exists, then the proportion of time that the process spends in state S _j is represented as:

${\pi }_j=\underset{t\to \top }{\mathsf{\text{lim}}}Pr\left(X\left(t\right)=S_j|X\left(0\right)=S_0\right)$

(13)

Therefore, we can define:

$\prod =\underset{t\to \top }{\mathsf{\text{lim}}}P\left(t\right)=\left[\begin{array}{cccc}\hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \\ \hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \end{array}\right]$

(14)

where N is the total number of states, subject to the normalisation condition $\sum _i{\pi }_i=1$. Then Eqn.12 becomes:

$\begin{array}{c}{\text{Π}}^{\prime }=\text{Π}*Q=\left[\begin{array}{cccc}\hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \\ \hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \end{array}\right]\left[\begin{array}{cccc}\hfill q_{S_1,S_1}\hfill & \hfill q_{S_1,S_2}\hfill & \hfill \cdots \hfill & \hfill q_{S_1,S_N}\hfill \\ \hfill q_{S_2,S_1}\hfill & \hfill q_{S_2,S_2}\hfill & \hfill \cdots \hfill & \hfill q_{S_2,S_N}\hfill \\ \hfill ⋮\hfill & \hfill ⋮\hfill & \hfill ⋮\hfill \\ \hfill q_{S_N,S_1}\hfill & \hfill q_{S_N,S_2}\hfill & \hfill \cdots \hfill & \hfill q_{S_N,S_N}\hfill \end{array}\right]=0\\ \left[\begin{array}{cccc}\hfill {\pi }_1\hfill & \hfill {\pi }_2\hfill & \hfill \cdots \hfill & \hfill {\pi }_N\hfill \end{array}\right]\left[\begin{array}{c}\hfill q_{S_1,S_i}\hfill \\ \hfill q_{S_2,S_i}\hfill \\ \hfill ⋮\hfill \\ \hfill q_{S_N,S_i}\hfill \end{array}\right]=0\mathsf{\text{ for }}i=1,\dots ,N\end{array}$

(15)

As $q_{S_i,S_i}=-\sum _{j\ne i}{q}_{S_i,S_j}$, then Equ.15 can be simplified to:

$\underset{\mathsf{\text{flux}}\mathsf{\text{out}}\mathsf{\text{of}}\mathsf{\text{state}}{S}_i}{\underset{}{{\pi }_i*\sum _{j\ne i}{q}_{S_i,S_j}}}=\underset{\mathsf{\text{flux}}\mathsf{\text{into}}\mathsf{\text{state}}{S}_{\mathsf{\text{i}}}}{\underset{}{\sum _{j\ne i}{\pi }_j*q_{S_j,S_i}}}$

(16)

We can therefore use Eqn.10,16 to obtain steady state distribution.

Here is a simple example to explain the above process. Suppose the system behaviour can be modeled by a two state Markov process whose state transition diagram is shown in Figure 3. Then the generator matrix is as follows:

$Q=\left[\begin{array}{cc}\hfill -7\hfill & \hfill 7\hfill \\ \hfill 6\hfill & \hfill -6\hfill \end{array}\right]$

(17)

We can therefore get steady state distribution from:

$\begin{array}{c}{\pi }_1*7={\pi }_2*6\\ {\pi }_1+{\pi }_2=1\end{array}$

(18)

which is:

$\pi =\left[{\pi }_1,{\pi }_2\right]=\left[0.46,0.54\right]$

(19)

Calculating the resource utilisation and system throughput

The steady state distribution π _i stands for the proportion of time that the process spends in state S _i . By finding in which state S _i the resource is busy, we can calculate the resource utilisation. For example, if the resource at place _i is busy in state S₁, S₃ and S₄, then its utilisation can be obtained by π(Resource_busy _placei ) = π(1) + π(3) + π(4). If this resource has multiple copies, we can calculate its average utilisation by

$Utilisatio{n}_{Resourc{e}_{placei}}=\left(\sum _{j=1}^{Nu{m}_{Resourcei}}\pi {\left(Resource\text{\_}bus{y}_{placei}\right)}_j\right)/Nu{m}_{Resourc{e}_{placei}}$

(20)

where $Nu{m}_{Resourc{e}_{placei}}$ is the total number of Resource _placei copies.

In the clinical pathway model, suppose that there is a patient who completes the activity with the action type λ leaves the hospital. The throughput of pathway is estimated by the expected number of completed activities with action type λ. Given the definition of resource Resource _placeN in the last place of the care process as follows:

$\begin{array}{cc}\hfill Resourc{e}_{-}idl{e}_{placeN}& \stackrel{def}{=}\left(\sigma ,rat{e}_{\sigma }\right).Resourc{e}_{-}bus{y}_{placeN}\hfill \\ \hfill Resourc{e}_{-}bus{y}_{placeN}& \stackrel{def}{=}\left(\lambda ,rat{e}_{\lambda }\right).Resourc{e}_{-}idl{e}_{placeN}\hfill \end{array}$

(21)

where the completion of activity (λ, rate_λ) enables the resource to become idle. Then the throughput T can be represented by:

$T=rat{e}_{\lambda }*Utilisatio{n}_{Resource\text{\_}bus{y}_{placeN}}$

(22)

showing that the system throughput is associated with the utilisation of the resource whose state transition depends on the type λ activity.

Let us introduce the CPP model for a simple clinical pathway (see Figure 4) to show how to analyze the resource utilisation and throughput. Suppose the CPP model is defined as follows:

(23)

The equilibrium state can only be maintained if the Markov process is irreducible that every state can be reached from all other states. Therefore CPP constructs a cyclic model in which the discharged patient will return to its initial state Patient_{place 0}. The state space of this simple CPP model is shown in Figure 5, where each state is represented by a five-tuple. For example, in state S₅: (0, 1, i, b, i), the first element 0 represents the first patient is in the place₀; the second element 1 represents the second patient is in the place₁; the third element i represents the first waiting room is idle; the fourth element b represents the second waiting room is busy; and the last element i represents the resource in the place₁ is idle.

The generator matix Q has the following form:

$Q=\left[\begin{array}{cccccccccc}\hfill -2\hfill & \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0.5\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1001\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1001\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 1000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1001\hfill & \hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1000\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1001\hfill & \hfill 1\hfill & \hfill 0\hfill & \hfill 1000\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -1000\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 500\hfill & \hfill 500\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -4\hfill & \hfill 0\hfill & \hfill 0.5\hfill & \hfill 0.5\hfill \\ \hfill 3\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -4\hfill & \hfill 0.5\hfill & \hfill 0.5\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 1.5\hfill & \hfill 0\hfill & \hfill 1.5\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -3\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 1.5\hfill & \hfill 0\hfill & \hfill 1.5\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -3\hfill \end{array}\right]$

(24)

from which we can get the steady state distribution:

$\begin{array}{c}\left[{\pi }_1,{\pi }_2,{\pi }_3,{\pi }_4,{\pi }_5,{\pi }_6,{\pi }_7,{\pi }_8,{\pi }_9,{\pi }_{10}\right]\hfill \\ =\left[0.5284,3.523*10^{-4},3.523*10^{-4},3.523*10^{-4},3.523*10^{-4},1.409*10^{-6},0.1761,0.1761,0.05895,0.05895\right]\hfill \end{array}$

(25)

Therefore, we can obtain resource utilisation and throughput:

$\begin{array}{cc}\hfill Utilisatio{n}_{Resourc{e}_{place1}}& =\pi \left(Resourc{e}_{-}bus{y}_{place1}\right)={\pi }_7+{\pi }_8+{\pi }_9+{\pi }_{10}=0.47\hfill \\ \hfill T& =3*Utilisatio{n}_{Resourc{e}_{place1}}=1.41\hfill \end{array}$

(26)

We can use PEPA eclipse plugin [30] to simulate the developed CPP model, from which resource utilization and throughput can be directly obtained.

Critical assessments-state explosion problem

When there are many patients and multiple copies of resources involved in the system description part, the state explosion problem occurs. For instance, when there are eight patients, eight waiting rooms and three resources in Figure 4 clinical pathway, the number of states increases from 10 to 75582. Then the dimension of generator matrix Q becomes 75582 × 75582. Hence, calculating steady state distribution by solving Eqn.16 turns out to be computational intensive and requires large storage space. To address these problems, this paper uses two methods which are state aggregation and fluid analysis.

It has been pointed out in [31, 32], by exploiting the strong equivalence relations, we can generate the aggregated CTMC where the number of states can be reduced. Consider again the simple CPP model in Eqn.23. As the states S₂, S₃, S₄ and S₅ shows that there is one patient in place 1, one patient in place 0, one waiting room is busy and the resource is idle, the aggregated state $S_2^{\prime }$ can be used to represent them without considering which patient is in place 1 and which waiting room is occupied:

(27)

Similarly, we use $S_4^{{}^{\prime }}$ to represent S₇ and S₈, and $S_5^{{}^{\prime }}$ to represent S₉ and S₁₀:

(28)

Therefore, the number of aggregated states is reduced to 5. The state transition diagram is shown in Figure 6, and the generator matrix Q' is as follows:

$Q^{\prime }=\left[\begin{array}{ccccc}\hfill -2\hfill & \hfill 2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill -1001\hfill & \hfill 1\hfill & \hfill 1000\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -1000\hfill & \hfill 1000\hfill & \hfill 0\hfill \\ \hfill 3\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -4\hfill & \hfill 1\hfill \\ \hfill 0\hfill & \hfill 3\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -3\hfill \end{array}\right]$

(29)

We can get steady state distribution as:

$\pi =\left[{\pi }_1,{\pi }_2,{\pi }_3,{\pi }_4,{\pi }_5\right]=\left[0.5284,0.001409,1.409*10^{-6},0.3523,0.1179\right]$

(30)

and further obtain resource utilisation and throughput:

$\begin{array}{cc}\hfill Utilisatio{n}_{Resourc{e}_{place1}}& =\pi \left(Resource\text{\_}bus{y}_{place1}\right)={\pi }_4+{\pi }_5=0.47\hfill \\ \hfill T& =3*Utilisatio{n}_{Resourc{e}_{place1}}=1.41\hfill \end{array}$

(31)

which are similar with Eqn.26.

Although the number of states can be reduced by using the aggregation method (implemented in PEPA eclipse plugin), when the amounts of patients and resources become extremely large, the state aggregation method cannot solve the state explosion problem well. Hence, fluid analysis of stochastic process model is used in [27] to analyze systems of size 10¹⁰⁰⁰ states and beyond. This approach approximates the state space by using a set of ordinary differential equations to describe the time evolution of state components. Assume that, the mean number of component P _i at time t is denoted as N (P _i (t)), then a differential equation can be used to represent the changes in N(P _i (t)) as:

$N{\left(P_i\left(t\right)\right)}^{\prime }=-\sum _{j:P_i\underrightarrow{(\alpha .,)}{P}_j}\mathsf{\text{rate of }}\alpha \mathsf{\text{ - action}}\mathsf{\text{leaving }}{P}_i+\sum _{j:P_j\underrightarrow{(\beta .,)}{P}_i}\mathsf{\text{rate of }}\beta \mathsf{\text{ - action}}\mathsf{\text{entering }}{P}_i$

(32)

A set of differential equations to present the CPP model in Eqn.23 is as follows:

$\begin{array}{cc}\hfill N{\left(Patien{t}_{place0}\left(t\right)\right)}^{\prime }& =-1*min\left(N\left(Patien{t}_{place0}\left(t\right)\right),N\left(Wait\text{\_}room0\text{\_}idle\left(t\right)\right)\right)\hfill \\ +3*N\left(Resource\text{\_}bus{y}_{place1}\left(t\right)\right)\hfill \\ \hfill N{\left(Patien{t}_{place1}\left(t\right)\right)}^{\prime }& =1*min\left(N\left(Patien{t}_{place0}\left(t\right)\right),N\left(Wait\text{\_}room0\text{\_}idle\left(t\right)\right)\right)\hfill \\ -3*N\left(Resource\text{\_}bus{y}_{place1}\left(t\right)\right)\hfill \\ \hfill N{\left(Wait\text{\_}room0\text{\_}idle\left(t\right)\right))}^{\prime }& =-1*min\left(N\left(Wait\text{\_}room0\text{\_}idle\left(t\right)\right)\right),N\left(Patien{t}_{place0}\left(t\right)\right))\hfill \\ +1000*N\left(Resourc{e}_{idle}\left(t\right)\right),N\left(\left(Wait\text{\_}room\text{\_}busy\left(t\right)\right)\right)\hfill \\ \hfill N{\left(Wait\text{\_}room0\text{\_}busy\left(t\right)\right))}^{\prime }& =1*min\left(N\left(Wait\text{\_}room0\text{\_}idle\left(t\right)\right)\right),N\left(Patien{t}_{place0}\left(t\right)\right))\hfill \\ -1000*N\left(Resourc{e}_{idle}\left(t\right)\right),N\left(\left(Wait\text{\_}room0\text{\_}busy\left(t\right)\right)\right)\hfill \\ \hfill N{\left(Resource\text{\_}idl{e}_{place1}\left(t\right)\right)}^{\prime }& =3*N\left(Resource\text{\_}bus{y}_{place1}\left(t\right)\right)\hfill \\ -1000*min\left(N\left(Resource\text{\_}idl{e}_{place1}\left(t\right)\right),N\left(Wait\text{\_}room0\text{\_}busy\left(t\right)\right)\right)\hfill \\ \hfill N{\left(Resource\text{\_}bus{y}_{place1}\left(t\right)\right)}^{\prime }& =-3*N\left(Resource\text{\_}bus{y}_{place1}\left(t\right)\right)\hfill \\ +1000*min\left(N\left(Resource\text{\_}idl{e}_{place1}\left(t\right)\right),N\left(Wait\text{\_}room0\text{\_}busy\left(t\right)\right)\right)\hfill \end{array}$

(33)

where the initial conditions are:

$\begin{array}{cc}\hfill N\mathsf{\text{(}}Patien{t}_{place0}\left(0\right))& \hfill \mathsf{\text{ = }}\mathsf{\text{2}}\\ \hfill N\mathsf{\text{(}}Patien{t}_{place1}\left(0\right))& \hfill \mathsf{\text{ = }}\mathsf{\text{0}}\\ \hfill N\left(Wai{t}_{-}room{0}_{-}idle\left(0\right)\right)& \hfill =2\\ \hfill N\mathsf{\text{(}}Wai{t}_{-}room{\mathsf{\text{0}}}_{-}busy(0))& \hfill =\mathsf{\text{0}}\\ \hfill N\mathsf{\text{(}}Resourc{e}_{-}idl{e}_{place\mathsf{\text{1}}}\left(0\right))& \hfill =\mathsf{\text{1}}\\ \hfill N\mathsf{\text{(}}Resourc{e}_{-}bus{y}_{place\mathsf{\text{1}}}(0))& \hfill =\mathsf{\text{0}}\end{array}$

(34)

By solving these equations, we can obtain the time evolution of each component. The details of fluid analysis can be found in [33]. By exploring the time evolution of P _i , we can find the mean time that P _i reaches its maximum. For example, Figure 7 shows N(Resource_busy_{place 1}(t)), whose maximum can be reached within two days. With the introduction of fluid analysis, state explosion problem can be addressed.

Calculating passage time

Based on fluid analysis, the mean passage time of a patient to be discharged can be estimated by using a stochastic probe [34]. For example, if we are interested in the time by which a Patient_{place 0}component has done its first α activity, we can attach a probe to Patient_{place 0}to remember whether it has performed α activity. The probe can be defined as follows:

$\begin{array}{cc}\hfill NotFinished& \stackrel{def}{=}\left(\alpha ,3\right).Finished\hfill \\ \hfill Finished& \stackrel{def}{=}\left(\alpha ,3\right).Finished\hfill \end{array}$

(35)

We can replace the Patient_{place 0}component in the system description part of the CPP model by the synchronised component . The modified differential equations are:

(36)

where the initial conditions are:

(37)

By summing the counts of all patient components of the form , we can get an approximation to the cumulative density function (CDF) for the time it takes for an individual Patient_{place 0}component to perform its first α action (see Figure 8). In this paper, we use a tool, Grouped PEPA Analyser, developed in [35] to simulate the CPP model. This tool can estimate resource utilisation and passage time based on fluid analysis.

Experiment setup

This paper applies CPP to model the main stream of the pathway where only the patients who can finally go to the HASU are analyzed. Hence the pathway in Figure 1 is simplified as it is shown in Figure 9.

The CPP model of the stroke clinical pathway

In order to get the average patient incoming rate to 4, the parameter Num _patient is set to 1000 and r _income is 0.004.

Three performance analysis scenarios

Three scenarios are tested in this subsection to show the application of the CPP model in performance analysis. The CPP model can detect the optimal resource allocation, estimate the passage time and determine the maximum throughput of the healthcare system.

Scenario 1: Resource optimization

Initially there are 3 stroke teams, 3 scanners and 20 beds available in the stroke clinical pathway. As the average patient input is around 4, it may happen that the resource utilisation is low. Hence, we investigate the influence of resource quantity on system throughput. By varying the numbers of stroke teams, scanners and beds, we can find from Figure 10 that the largest throughput that can be achieved is 3.95 patients/day, slightly smaller than 4 patients/day. We cannot have the maximum throughput value exactly equal to the number of input patients. It can be explained from the definition of parallel input patients Patient _Home [Num _patient ], where Num _patient is 1000 and the input rate is 0.004. Theoretically, the number of input patients per day should be 1000 * 0.004 = 4. On the first day, 4 patients on average fall ill and are treated sequentially through the clinical pathway. As the mean time for the patient staying in HASU department is 2 days, it is possible that these 4 patients are still in the pathway on the next day. Hence, they cannot return to the state Patient _Home on the second day, meaning that the number of patients at Patient _Home is slightly smaller than 1000. Therefore, there is a trivial difference between the maximum throughput and the number of input patients. However, this small difference does not influence the performance analysis process.

Figure 10 is used to show how our performance analysis result guides medical staff to find optimal resource allocation, i.e., the amount of resource is kept minimal while still maintaining the maximum throughput in the clinical pathway. The three 3D graphs look similar and they actually demonstrate the process of reducing redundant resources while keeping the throughput unchanged. In these 3D graphs, we give two examples. In the first example, the numbers of beds and scanners are reduced from 20/3 to 8/1, respectively. This example shows that twelve redundant bed resources and two scanner resources can be removed and the clinical pathway's throughput is not influenced. The second example further reduces two redundant stroke team resources and keeps the throughput unchanged.

We can obtain the optimal number of each resource in Table 3 and find that the resource utilisation can be increased while the system throughput is maintained. It is notable that the CPP model can find the optimal resource allocation with varying number of input patients. For example, with an input of 4 patients per day, the optimal resource allocation is 1 stroke team, 1 scanner and 8 beds; if the input rate increases to 0.006, meaning that there are 6 patients falling ill every day, the optimal resource allocation is found to be 1 stroke team, 1 scanner and 12 beds.

	Stroke team		Scanner		Bed
	Number	Utilisation	Number	Utilisation	Number	Utilisation
Initial parameter set	3	0.038	3	0.153	20	0.368
Optimal parameter set	1	0.115	1	0.46	8	0.92

Scenario 2: Passage time estimation

Besides estimating the system throughput and resource utilisation, CPP model can further predict the mean passage time of a patient with the introduction of a stochastic probe as follows:

$\begin{array}{cc}\hfill NotFinished& \stackrel{def}{=}\left(treat\text{\_}HASU,r_{treat\text{\_}HASU}\right).Finished;\hfill \\ \hfill Finished& \stackrel{def}{=}\left(treat\text{\_}HASU,r_{treat\text{\_}HASU}\right).Finished;\hfill \end{array}$

(41)

The state NotFinished becomes Finished once the activity whose action type is treat_HASU is completed. When we attach this probe to the patient state component as follows:

(42)

the mean time of a patient from falling ill to leaving the HASU can be obtained. Figure 11 shows the probability that a patient is discharged from the HASU within a specific time duration when there is 1 stroke team, 1 scanner and 8 beds available. We can see that 90% of patients can be discharged within 2.5 days and only 14% of patients can be discharged within 1 day. Estimation of the whole passage time not only tells patients when they can be discharged, but also informs doctors the resource occupation time.

Scenario 3: Maximum input estimation

We can also examine the maximum throughput with different parameter settings. For example, if the pathway model uses the initial parameters, in Figure 12 we can see that when the number of input patients increases from 1 to 100, the maximum throughput is around 10 patients/day. When the input number is 10, the throughput reaches 9.74 patients/day, slightly smaller than the maximum value (this difference has been discussed in the Resource optimization subsection). Any further increase in the number of input patients has trivial contribution to the system throughput.

Therefore, by simulating the CPP model, we can find the maximum number of input patients that can be supported by the healthcare system. This maximum input estimation is significant for hospital to determine whether it can accept more stroke patients or not. For example, suppose there are already 10 patients on average coming to the hospital from the surrounding area. If the national health community asks whether this hospital can serve patients from larger area, meaning that more than 10 patients will arrive every day, by estimating the maximum input this hospital can determine whether this is possible and whether more resources are required to support the increased number of patients.