**Table 9.** Analysis based on k-fold.

**Table 10.** Comparative analysis.

Using cluster size = 9, and considering the Cleveland database, the proposed Taylor-BSA–DBN showed maximal accuracy of 0.934, which is 16.92%, 11.13%, and 3.96%, better than the existing methods, such as SVM, NB, and DBN, respectively. In the existing methods, the DBN offers maximum sensitivity of 0.913, but the proposed method is 3.89% better than the existing DBN. The proposed method has a maximum specificity of 0.903. The percentage of improvement of the proposed method with the existing methods, such as SVM, NB, and DBN, is 23.15%, 15.28%, and 3.10%, respectively. Considering the Hungarian database, the proposed Taylor-BSA–DBN showed maximal accuracy of 0.902, maximal sensitivity of 0.909, and maximal specificity of 0.893. Considering the Switzerland database, the proposed Taylor-BSA–DBN showed maximal accuracy of 0.840, which is 19.17%, 10.12%, and 2.38%, better than the existing methods, such as SVM, NB, and DBN, respectively. Similarly, the proposed system has a maximum sensitivity of 0.846. The percentage of improvement of the proposed system sensitivity with the existing methods, such as SVM, NB, and DBN is 19.74%, 11.35%, and 1.89%, respectively. Likewise, the proposed Taylor-BSA–DBN showed maximal specificity of 0.833.

Table 11 shows the computational time of the proposed system and the existing methods, such as SVM, NB, and DBN, in which the proposed Taylor-BSA–DBN has a minimum computation time of 6.31 sec.

Table 12 shows the statistical analysis of the proposed work and the existing methods based on mean and variance.

**Table 12.** Statistical Analysis.

#### **5. Conclusions**

Contemporary medicine depends on a huge amount of information contained in medical databases. The obtainability of large medical data leads to the requirement of effective data analysis tools for extracting constructive knowledge. This paper proposes a novel, fully automated DBN for heart disease diagnosis using medical data. The proposed Taylor-BSA is employed to train DBN. The proposed Taylor-BSA is designed by combining the Taylor series and BSA algorithm, which can be utilized for finding the optimal weights for establishing effective medical data classification. Here, the sparse-FCM is employed for selecting significant features. The incorporation of sparse FCM for the feature selection process provides more benefits for interpreting the models, as this sparse technique provides important features for detection, and can be utilized for handling high dimensional data. The obtained selected features are fed to DBN, which is trained by the proposed Taylor-BSA. The proposed Taylor-BSA is designed by integrating the Taylor series and BSA in order to generate optimal weights for classification. The proposed Taylor-BSA–DBN outperformed other methods with maximal accuracy of 93.4%, maximal sensitivity of 95%, and maximal specificity of 90.3%, respectively. The proposed method does not classify the type of heart disease. In the future, other medical data classification datasets will be employed for computing efficiency of the proposed method. In addition, the proposed system will be further improved to classify heart diseases, such ascongenital heart disease, coronary artery disease, and arrhythmia.

**Author Contributions:** Conceptualization, A.M.A. methodology, A.M.A.; software, A.M.A.; validation, A.M.A.; resources, A.M.A.; data curation, A.M.A.; writing—original draft preparation, A.M.A. and W.M.N.W.Z.; writing—review and editing, A.M.A. and W.M.N.W.Z.; visualization, A.M.A. and W.M.N.W.Z.; supervision, W.M.N.W.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Decision Support System for Elective Surgery Scheduling under Uncertain Durations**

**Daniel Clavel 1, Cristian Mahulea 1,\*, Jorge Albareda 2 and Manuel Silva 1**

Received: 22 January 2020; Accepted: 9 March 2020; Published: 12 March 2020

**Featured Application: A software tool for helping doctors in the automatic scheduling of elective patients is proposed. The tool has been developed considering the organizational structure of Spanish hospitals. However, it may be easily extended to any hospital (public or private) working under operation room block booking criterion.**

**Abstract:** The *operation room* (OR) is one of the most expensive material resources in hospitals. Additionally, the demand for surgical service is increasing due to the aging population, while the number of surgical interventions performed is stagnated because of budget reasons. In this context, the importance of improving the efficiency of the surgical service is accentuated. The main objective of this work is to propose and to evaluate a *Decision Support System* (DSS) for helping medical staff in the automatic scheduling of elective patients, improving the efficiency of medical teams' work. First, the scheduling criteria are fixed and then the scheduling problem of elective patients is approached by a mathematical programming model. A heuristic algorithm is proposed and included in the DSS. Moreover, other different features are implemented in a software tool with a friendly user interface, called CIPLAN. Considering realistic data, a simulation comparison of the scheduling obtained using the approach presented in this paper and other similar approaches in the bibliography is shown and analyzed. On the other hand, a case study considering real data provided by the Orthopedic Surgical Department (OSD) of the "Lozano Blesa" hospital in Zaragoza (HCU) is proposed. The simulation results show that the approach presented here obtains similar occupation rates and similar confidence levels of not exceeding the available time than approaches in the bibliography. However, from the point of view of respecting the order of the patients in the waiting list, the approach in this paper obtains scheduling much more ordered. In the case of the Orthopedic Surgical Department of the "Lozano Blesa" hospital in Zaragoza, the occupation rate may be increased by 2.83%, which represents a saving of 110,000 euros per year. Moreover, medical doctors (who use this tool) consider CIPLAN as an intuitive, rapid and efficient software solution that can make easier the corresponding task.

**Keywords:** decision support system; heuristic approach; surgery scheduling; software tool; case study

#### **1. Introduction**

The *Operation Room* (OR) is one of the most expensive material resources of the hospitals (in Spain being about 40% of the total hospitalization cost [1]) and approximately 60% of patients need it at some point during their hospital stay [2]. On the other hand, the demand for surgical services is increasing due to the aging population. In Spain, the average life expectancy is about 81.8 years being one of the highest in Europe [3].

Moreover, as a result of the economic crisis in 2008, the number of surgical interventions performed annually in Spanish public hospitals has stagnated at 3.5 million. This situation has also contributed to a continuous increment in the number of patients in the surgical waiting lists, reaching 600,000 in 2017 [4]. In some cases, these huge lists lead into waiting times greater than six months, which is becoming a serious problem for the healthcare system. Being difficult to increase the surgical resources, the importance of improving the efficiency is accentuated. New management techniques, scheduling methods, and specific information systems could help to improve the efficiency of the surgical services and, consequently, could reduce the surgical waiting lists.

The main objective of this work is to provide a DSS for managing surgical services in hospitals. In particular, a software tool (called CIPLAN from the Spanish words *PLANificación de CIrugías*) is presented that may be used by medical doctors to automatically perform a scheduling of the patients. It is assumed that the uncertain parameters (surgery duration and cleaning time) follow a normal distribution and, consequently, the expected total duration of a surgical block also follows a normal distribution. In this way, the probability of exceeding the working time can be introduced as a capacity chance constraint. Even if the lognormal distribution may give better approximations of the surgery times [5], the optimization problem has nonlinear constraints, being difficult to solve. Since, in general, acceptable error is obtained for normal distribution, this work assumes this type of distribution.

The following three criteria have been considered for the scheduling problem.

*C1. Maximize the occupation of the ORs.* Since the ORs with their associated cost constitute the bottleneck of the system (they can be open only during a given period of time every working day and their cost is high), the first scheduling criterion is to maximize the occupation of the ORs. In this way, the underutilization of the ORs is prevented and the use of resources is maximized.

*C2. Ensure with a minimum confidence level that the total available time in a time block is not going to be exceeded.* On contrary to *C1*, an over-scheduling of the ORs is not desired since the last scheduled surgeries can be canceled if the medical personnel believes that the time of the time block will be exceeded. Among other problems, such cancellations reduce the *real* occupation rate in the ORs.

*C3. Respect as much as possible the order of the patients in the waiting list.* Due to ethical and common sense reasons, patients with the same urgency level should spend a similar time in the waiting list. For that, once the patients in a waiting list are ordered depending on the urgency level and the admission date, the third scheduling criterion is to respect as much as possible the order of the patients in the waiting list. For example, in the Orthopedic Department of the HCU, each team has assigned two time blocks per week for the elective patients, so the difference of scheduling a patient in the first time block with the sixth one is 21 days.

*Mathematical programming* is an optimization technique [6] commonly used in managing hospitals [7], not only for surgery scheduling problems [8], but also for other problems such as bed management [9], nurse scheduling [10], or management policy decisions [11]. This paper proposes first an *Mixed Integer Quadratic Constrained Programming* problem for performing the scheduling of patients based on the three criteria presented before. It is well known that the computational complexity to obtain the optimal solution is high and a heuristic method is proposed called *Specific Heuristic Algorithm* (SHA) [12]. The software tool, in addition to helping doctors in the scheduling task, has also been developed for managing medical doctors, medical teams, OR timetable, and patients.

A preliminary version of this work has been reported in [12] where a first version of the SHA has been proposed. In this paper, the heuristic SHA is improved by adding a post arranging of the time scheduled blocks. This improvement provides solutions respecting the order of the patients in the waiting list more. Additionally, the DSS is explained (not only the part related to surgery scheduling but also related to management tasks), and the solution is compared using real data from the HCU with two other methods in literature.

#### **2. Related Work**

In literature, many researchers studied the problem of planning and scheduling of elective patients (see [8,13] and the references herein for an overview while for an overview of the problems see [14]). Almost all the works consider criterion *C1* presented before. Criterion *C2* (of using chance constraints for the scheduling) is considered also in some contributions presented in the last years. For example, in [15], integer stochastic formulations are proposed for two sub-problems in which the original problem is decomposed. The approach to solve them is based on neighborhood search techniques and Monte Carlo simulation. In [16], a two-stage robust optimization model is proposed that is solved based on a column-and-constraint generation method. A similar method is used by the authors in [17] to deal with uncertain surgery duration and emergency demand, while, in [18], the authors use particle swarm optimization for the scheduling of surgeries.

The approaches presented in [19,20] are the most related to the one presented in this paper. The authors in [19] present an optimization framework for batch scheduling within a block booking system that maximizes the expected utilization of ORs subject to a set of probabilistic capacity constraints. They propose an algorithm that iteratively solves a series of mixed-integer programs that are based on a normal approximation of cumulative surgery duration. In [20], constructive and local search heuristics for maximization of ORs utilization and minimization of the overcoming risk are proposed. In their model, to address the randomness of surgery processing times, a planned time slack is reserved in each scheduling block, which is a function of total mean and variance of surgeries assigned to the corresponding block. When determining an appropriate size of the planned slacks, the authors assume that the sum of surgery duration follows a normal distribution.

Both approaches in [19,20] require setting in advance the patients that are going to be scheduled in the next blocks. Therefore, in these approaches, all considered patients must be scheduled in one of the available blocks. For this, both approaches start with an initial scheduling obtained through the scheduling rule: *first-fit probabilistic*. Following this rule, sequentially each surgery is assigned to the first available block for which the probabilistic capacity constraint is satisfied after the assignment. Once the initial scheduling is obtained, the expected occupation rate of the first blocks is improved by rescheduling the surgeries. In this way, the last blocks are totally or partially released.

Unlike [19,20], the approach presented in this paper does not require to know in advance the set of patients that should be scheduled in the next surgical blocks, since any patient on the waiting list may or may not be scheduled. Another big difference with [19,20] consists of the usage of two balanced terms in the cost function that favor patients to be scheduled in an orderly manner at the same time that the expected occupation rate of the OR is maximized. An exhaustive numerical comparison using realistic data is performed and analyzed showing that the approach presented in this manuscript provides better solutions to the considered problem than the one in [19,20]. In particular, even if the results show similar occupation rates and similar confidence levels, the approach presented in this paper respect the preference order of the patients in the waiting list more. Therefore, criterion *C3* is better considered. In the case of the hospital in which the approach is implemented, the number of available time blocks per week for each medical team is equal to two, while, on average, three patients are scheduled in a time block, making it very important to keep, if possible, the order of patients in the waiting lists. This is one of the main differences of the approach presented in this paper with the other approaches.

DSS are used for managing tasks in different fields of applications, as for example in stormwater management [21], dental manufacturing production line [22], offshore wind farms [23] and healthcare systems [24,25]. Moreover, a DSS for OR scheduling was proposed in [26] based on a Mixed Integer Linear Programming (MILP) problem [27,28] that takes into account the preference order of the patients. However, this approach does not consider the variability on the uncertain parameters (surgery duration, cleaning time, etc.) and consequently can not manage the risk of exceeding the total time. The DSS proposed in this work uses an heuristic algorithm allowing to impose a minimum confidence level of not exceeding the total time in a surgical block. In this way, doctors not only can manage the risk but also can know in advance the probability that each block exceeds the available time. The DSS proposed in this paper can also perform the management of medical teams and OR time table. Moreover, it is possible to perform iterative scheduling which is useful if, for example, after obtaining a first scheduling, some patients can not attend the day of their surgeries. Furthermore, after a surgery is performed, the DSS automatically updates and customizes the average and standard deviation of the surgery duration. Finally, in order to provide a safe environment to work and better security of the medical data, different access levels are provided to every user.

The second main contribution of this paper is the presentation of a DSS that can help medical doctors not only for scheduling of surgeries, but also it can be used for management tasks of patients, medical teams, hospital resources, etc. Finally, the scheduling approach presented is compared with real historical data from the studied hospital, this being a third contribution of this manuscript.

#### **3. Preliminaries and Problem Definition**

Let us assume that each surgical team of a hospital department has some surgical blocks booked for performing the surgeries of their elective patients list. Thus, the problem of scheduling can be divided into different independent sub-problems: one for each medical team. In this process, some patients are assigned from the waiting list to the next time blocks previously booked. Currently, in many Spanish hospitals, the scheduling is performed manually, so the schedulers (normally medical doctors) are guided by their own intuition and experience to assign the patients.

#### *3.1. Scheduling Problem Statement*

In this subsection, the **Scheduling problem** is introduced formally. First, the notation, terminology, and assumptions are fixed and then the scheduling criteria are transformed to cost function and constraints.

Let S = {*s*1,*s*2, ... ,*s*|S|} be the set of surgery types that can be performed in the medical department and let us assume that the duration *d*(*sk*) of each type of surgery *sk* ∈ S is a random variable with normal probability density function (pdf) *d*(*sk*) = *N*(*μd*(*sk* ), *σd*(*sk* )), where *μd*(*sk* ) is the average and *σd*(*sk* ) is the standard deviation.

On the other hand, let *Dt* be the delay with respect to the starting time of a block and let *Ct* be the cleaning time duration of the OR after a surgery is performed. It is assumed that *Dt* and *Ct* are also random variables with normal pdf, i.e., *Dt* = *N*(*μDt*, *σDt*) and *Ct* = *N*(*μCt*, *σCt*).

Furthermore, let us consider W = {*w*1, *w*2, ... , *w*|W|} an ordered list of patients such that, if *wj* ∈ W, *j* is the preference order of the patient *wj* in the waiting list.

In addition, let us assume an ordered set of time blocks B = {*b*1, ..., *b*|B|} to schedule, where *b*|B| is the block corresponding to the latest date. Each block *bi* ∈ B has a fixed duration denoted by *l*(*bi*).

For each time block *bi* ∈ B to schedule, there exists a binary decision vector *Si* ∈ {0, 1}1×|W| with a dimension equal to the number of the patients in the waiting list. If *Si*[*j*] = 1, then surgery of patient *wj* should be performed in working day *i* ≤ |B|.

Let *Sbi* be the set of surgeries scheduled in the time block *bi*; then, the total duration *Tbi* of the time block *bi* can be computed by Equation (1). That is the sum of: (i) the delay with respect to the starting time (i.e., *Dt*), (ii) the duration of the surgeries scheduled in the time block *bi*, and (iii) the cleaning time duration *Ct* between two consecutive surgeries. Notice that the cleaning time after the last surgery is not considered:

$$T\_{b\_i} = Dt + \left( |S\_{b\_i}| - 1 \right) \cdot \mathbb{C}t \sum\_{s\_k \in S\_{b\_i}} \left[ d(s\_k) \right]. \tag{1}$$

Because all these variables are assumed to follow a normal distribution, the total duration *Tbi* of the time block *bi* also follows a normal distribution *Tbi* = *N*(*μTbi* , *σTbi* ).

The goal of the scheduling problem is the assignment of the patients from the waiting list W to the set of time blocks B considering the scheduling criteria C1, C2 and C3. Let us define the following further notations:

then:

Criterion C1 is formalized by objective **O1** given by Equation (2). This equation maximizes the expected surgery time in each block (*μw* · *Si*) giving more importance to the first ones (|B| + 1 − *i*):

$$\max \sum\_{i=1}^{|\mathcal{S}|} \left[ \mu\_{\text{tw}} \cdot \mathcal{S}\_{i} \cdot \left( |\mathcal{B}| + 1 - i \right) \right]. \tag{2}$$

Criterion C2 is imposed by the chance probability constraint given by Equation (3). It ensures with a minimum confidence level *Cl* that the total duration *Tbi* in each time block *bi* ∈ B does not exceed the available time *l*(*bi*):

$$P[T\_{b\_i} \le l(b\_i)] \ge \mathcal{C}l, \forall i = 1, 2, \dots, |\mathcal{B}|. \tag{3}$$

Finally, criterion C3 is formalized by the objective **O2** given by Equation (4). Here, by minimizing the sum of the preference order of the patients scheduled in each time block (*Po* · *Si*), preference to the first patients in the waiting list is given. Moreover, multiplying by the term (|B| + 1 − *i*), first patients should be scheduled in the first time blocks:

$$\min \sum\_{i=1}^{|\mathcal{B}|} \left[ \mathbf{P} \mathbf{o} \cdot \mathbf{S}\_i \cdot \left( |\mathcal{B}| + 1 - i \right) \right],\tag{4}$$

#### *3.2. Optimal Mathematical Programming Model*

A Mixed-Integer Quadratic Constrained Programming (MIQCP) model is proposed (5) for solving the patients **Scheduling problem**. In this model, the objectives given in Equations (2) and (4) are mixed in a linear cost function with two terms balanced by the value of the parameter *β*:

$$\begin{aligned} \min \quad & \sum\_{i=1}^{|\mathcal{B}|} \left[ -\mu\_w \cdot \mathbf{S}\_i \cdot (|\mathcal{B}| - i + 1) \right. \\ & \quad \pounds \cdot \mathbf{P}\_o \cdot \mathbf{S}\_i \cdot (|\mathcal{B}| - i + 1) \right] \end{aligned} $$

Subject to:

$$\begin{cases} P[T\_{b\_i} \le l(b\_i)] & \ge \quad \mathcal{C}l, \quad \forall i = 1, 2, \dots, |\mathcal{B}|. \\ \sum\_{i=1}^{|\mathcal{B}|} \mathcal{S}\_i[j] & \le \quad 1, \quad \forall j = 1, 2, \dots, |\mathcal{W}| \\ \mathcal{S}\_i \in \{0, 1\}^{1 \times |\mathcal{W}|}, & \forall i = 1, 2, \dots, |\mathcal{B}|. \end{cases}$$

The MIQCP model has two sets of constraints. The first one is related to the chance probability constraint given in Equation (3), while the second one imposes that each patient in the waiting list is scheduled at most once. Development of the chance constraint (3) quadratic constraints is obtained [29].

(5)

Moreover, the decision variables *Si* are binary, so the computational complexity of the MIQCP model is really high. In order to obtain scheduling in a reasonable time, heuristic approaches are needed.

#### **4. Proposed Solution**

#### *4.1. Specific Heuristic Algorithm*

Two heuristic approaches have been discussed for implementation in the DSS. The first one is a *Receding Horizon Strategy* (RHS) which obtains a sub-optimal solution sequentially (similar to [30]), while the second one is a *Specific Heuristic Algorithm* (SHA) based on list scheduling techniques [31–33].

The SHA for the scheduling of patients problem is composed of three parts: (i) a previous data analysis; (ii) the scheduling of the time blocks; and (iii) a re-assignment of time blocks.

Part (i): A previous data analysis. This part is composed of two steps executed sequentially once at the beginning of the scheduling.

*Step 1. Classify the set of surgeries* S in a given number *t* of disjoint subsets such that S = *S*1 ∪ ... ∪ *St*, *t* being an input parameter. This classification is done based on two conditions: (1) the average duration of all surgeries in a subset *Si* ∈ S, *i* = 1, 2, ... , *t* − 1 must be less than or equal to the average duration of all surgeries in subset *Si*+1 ∈ S, and (2) the expected number of patients in each subset *Si* ∈ S should be similar. In order to compute the expected number of patients for each *Si* ∈ S, the probability of occurrence of all surgeries in S should be known, and it is obtained using historical data of the studied department.

For example, if there exist three types of surgeries S = {*s*1,*s*2,*s*3} with *μd*(*s*1) < *μd*(*s*2) < *μd*(*s*3) and the occurrence probabilities are 0.25 for *s*1 and *s*2 and 0.5 for *s*3, a valid classification for *t* = 2 is S = {*S*1, *S*2} where *S*1 = {*s*1,*s*2} and *S*2 = {*s*3}. However, if *μd*(*s*1) < *μd*(*s*3) < *μd*(*s*2), a valid classification could be *S*1 = {*s*1,*s*3} and *S*2 = {*s*2} but also *S*1 = {*s*1} and *S*2 = {*s*2,*s*3} is a valid solution.

*Step 2: Obtain the set of possible scheduling types for each time block bi* ∈ B. In this step, all multisets from set S = {*S*1, *S*2, ... , *St*} are obtained from the possible scheduling type for the time block *bi*. Formally, a possible scheduling type is a tuple (S, *m*), where S is the underlying set of the multiset and *m* : S → N≥1 is a function giving the multiplicity of the elements in the multiset. A multiset belonging to the possible scheduling type of block *bi* should satisfy the chance constraint (3). In order to check that the chance constraint (3) is satisfied, the surgeries with lower average duration are considered.

For example, a possible scheduling type for time block *bi* could be {*Sk*, *Sk*, *Sj*} meaning that, in block *bi*, two surgeries from set *Sk* may be scheduled and another one from set *Sj*. If *s*1 *k* is the surgery type with the lowest average duration from *Sk* and *s*1 *j* is the surgery type with the lowest average duration from *Sj*, according to the total duration in (1), the minimum total duration of the time block *bi* is

$$T\_{b\_{\hat{\mathbf{y}}}} = \left(Dt + s\_k^1 + s\_k^1 + s\_{\hat{\jmath}}^1 + 2 \cdot Ct\right) ,$$

where *Dt* is the delay with respect to the starting time and *Ct* is the cleaning time duration after each surgery.

As observed before, *Tbi* is following a normal distribution with mean

$$
\mu\_{T\_{b\_{\bar{l}}}} = \mu\_{Dt} + 2 \cdot \mu\_d(s\_k^1) + \mu\_d(s\_{\bar{j}}^1) + 2 \cdot \mu\_{Ct},
$$

and standard deviation

$$
\sigma\_{T\_{\boldsymbol{b}\_{i}}} = \sqrt{\left(\sigma\_{Dt}\right)^{2} + 2\cdot \left(\sigma\_{d}(\mathbf{s}\_{k}^{1})\right)^{2} + \left(\sigma\_{d}(\mathbf{s}\_{\boldsymbol{\zeta}}^{1})\right)^{2} + 2\cdot \left(\sigma\_{\mathrm{Ct}}\right)^{2}}.
$$

Since {*Sk*, *Sk*, *Sj*} is a possible scheduling type, the chance constraint (3) is satisfied. In order to check this equation, based on normal distribution property [34], the following inequality

$$P[T\_{b\_i} \le l(b\_i)] \ge CI$$

is equivalent to

$$\frac{l(b\_i) - \mu\_{T\_{b\_i}}}{\sigma\_{T\_{b\_i}}} \ge V\_{Cl\prime}$$

where *VCl* is the value corresponding to the normal variable (*x* ∼ *N*(0, 1)) with a cumulative distribution function *Cl*, i.e., *P*[*x* ≤ *VCl*] = *Cl* and this value is tabulated.

Part (ii): Scheduling of the time blocks, consisting of the following four steps executed sequentially for each time block *bi* to schedule.

*Step 3: For each possible scheduling type of time block bi, obtain the sets of real scheduling*. Given a possible scheduling type, a real scheduling is composed by a set of patients with the types of surgeries equal to the possible scheduling types.

For example, considering the scheduling type {*Sk*, *Sk*, *Sj*}, a real scheduling is composed from two patients *w* and *w* on which a surgery *sk* ∈ *Sk* should be performed and one patient *w* on which a surgery *sj* ∈ *Sj* should be performed. The first real scheduling corresponding to {*Sk*, *Sk*, *Sj*} is chosen by taking the patients with a lower number of preference order (corresponding to the higher preference). However, all other combinations of patients satisfying the chosen scheduling type having no patients with preference order greater than the maximum preference order of the patient in the first real scheduling will be considered.

Let us assume, for example, the list of patients given in Table 1 with W = {*w*1, *w*2, *w*3, *w*4, *w*5, *w*6} such that patients *w*1, *w*3, and *w*5 have assigned the surgery *s*1, the patient *w*2 has assigned the surgery *s*2, the patient *w*4 the surgery *s*3, and the patient *w*6 the surgery *s*4. Assuming the set of surgeries S = {*S*1, *S*2, *S*3} with *S*1 = {*s*1,*s*2}, *S*2 = {*s*3} and *S*3 = {*s*4} and two possible scheduling types {*S*1, *S*1, *S*2} and {*S*1, *S*2, *S*3}, the first real scheduling considered for {*S*1, *S*1, *S*2} is {*w*1, *w*2, *w*4}. The maximum order of patients in this real scheduling is 4 (corresponding to *w*4), hence the following real scheduling will be considered as well: {*w*1, *w*3, *w*4} and {*w*2, *w*3, *w*4}. For the scheduling type {*S*1, *S*2, *S*3}, the first real scheduling is {*w*1, *w*4, *w*6}. Since the maximum preference order of patients is 6, the following real scheduling will be considered as well: {*w*2, *w*4, *w*6}, {*w*3, *w*4, *w*6} and {*w*4, *w*5, *w*6}. All of this real scheduling is shown in Table 2.

**Table 1.** Surgery data and waiting list of patients used for explaining part 2 of the SHA.

*Step 4: Evaluate the real scheduling and select the best one for each scheduling type*. First, the real scheduling not fulfilling the chance constraint (3) is removed and the check is similar to the one in Step 2. For the remaining real scheduling, the expected occupation rate *r* is computed by summing the average duration of all surgeries in the real scheduling divided by the time block duration *l*(*bi*) (multiplied by 100 if percentage). Furthermore, the average preference order *Ap* is computed by summing the preference order of all patients in the real scheduling divided by the number of patients in the real scheduling. The real scheduling with the minimum value of the following fitness function is chosen:

$$H = (Ap - Min\_{Ap}) \cdot \beta + (Max\_r - r),\tag{6}$$

where *MinAp* is the minimum average preference order of all real scheduling; *Maxr* is the maximum occupation rate *r* of all real scheduling; and *β* is a balancing parameter between the two terms (one for preference order and one for occupation rate).

For example, in Table 2, the computations of the occupation rates are given, preference order, and fitness function (fifth column) for the set of real scheduling considered before. In this example, a value of *β* = 2.6 is used and the time block duration are assumed to be *l*(*bi*) = 420 (min), corresponding to 7 h. In this step, the real scheduling *w*1*w*2*w*4 is chosen for the possible scheduling type *S*1*S*1*S*2 and *w*2*w*4*w*6 is chosen for *S*1*S*2*S*3.

*Step 5: Evaluate the best real scheduling chosen in Step 4 (one for each scheduling type) and select the best one for the current time block bi*. From the set of selected real scheduling in the previous step, the one with minimum value of the fitness function (6) is chosen. Notice that the value of the fitness function of a real scheduling in this step is, in general, different by the one computed in step 4. In particular, *MinAp* and *Maxr* have different values since are computed based only on the best real scheduling computed in step 4.

The sixth column of Table 2 is showing the computation of the fitness function on the two selected real scheduling instances on the previous step. Hence, in this step, *MinAp* = 2.33 corresponding to *w*1*w*2*w*4 while *Maxr* = 88.09 corresponding to *w*2*w*4*w*6. According to them, patients *w*2*w*4*w*6 are selected for the time block *l*(*bi*).

*Step 6: Remove the scheduled patients*. The patients scheduled in Step 5 are removed from the waiting list, and the procedure is iterated for time block *bi*+1 from Step 3.

By using the case study in the HCU hospital, it has been observed that values of *t* = 3 and of *β* = 2.6 are appropriate. However, for different medical departments, a previous study is necessary in order to better choose these design parameters.

Part (iii): Sorting the time blocks with the same available time. First, the time blocks are grouped depending on their available time *l*(*bi*). Then, the time blocks of each group are sorted in ascending order according to the average preference order of the scheduled patients (indicator *Ap*). At the end, for each date, a time block is assigned sequentially starting with the closest one.

#### *4.2. Decision Support System*

In this section, the main features of the proposed DSS are discussed. Each doctor in the hospital department has his own list of patients and the waiting list of a team is composed by the merged list of doctors belonging to the team. Each team must schedule the patients from its waiting list during the time blocks previously booked by the head of the department. Thus, the main objective of the DSS is to help medical doctors perform a fast efficient and dynamic scheduling of the elective patients. The DSS uses as a core for the scheduling of patients the SHA approach presented in the previous subsection. However, other features are also included which enable (a) managing the medical teams and the OR time-table; (b) updating the waiting lists; (c) dynamic planning; and (d) improving the problem parameter by updating the surgeries duration.

*Manage medical team and OR time-table.* Normally, the medical teams are composed of the same medical doctors. However, sometimes medical doctors could move from one team to another. The head of the department is responsible for updating the medical teams and including this information in the DSS. Moreover, the head of the department should define the ORs time-table in the DSS. The OR, the date, the starting time, and the ending time of each time block assigned to each team is the information that must be included in the DSS.

*Updating the waiting list.* Generally, the new arrived patients are included at the end of the waiting list. However, depending on medical criteria, the surgeon could decide to advance the patient into a higher position in the waiting list. The DSS automatically orders the waiting list of patients. There are two parameters influencing directly the position of the patient in the waiting list.

The final score is obtained as follows:

$$S\_T = p \cdot S\_1 + S\_2 \tag{7}$$

Finally, the patients are ordered according to their total score. The patient with the highest total score will be the first, while the patient with the lowest punctuation will be the last one.

*Iterative planning.* The scheduler of each medical team performs the scheduling for the next |B| time blocks (this is done by using the SHA). The secretary calls the patients scheduled and asks them for their availability on the scheduled date. The secretary gives back this information to the team scheduler who should schedule the empty gaps again. This process is repeated until the |B| time blocks are completely scheduled with all patients confirmed.

During the iterative scheduling, if a patient confirms the attendance, she/he is fixed in the corresponding time block. However, in the case that a patient cannot be contacted or she/he cannot be hospitalized, then the patient is not considered for the scheduling of the next |B| time blocks.

Thus, in the next iteration, the patients previously confirmed and the patients who cannot be hospitalized are not considered in the waiting list. The SHA schedules new patients in the gaps of the uncompleted time blocks.

*Updating and customizing the average duration.* The average duration and standard deviation of each type of surgery should be computed using historical data from the hospital department. Considering two years, a sufficiently high number of surgeries of each type is obtained, and the average duration is representative. However, depending on the skills and experience of the surgeons, significant differences could exist between their average duration. Moreover, for the same medical doctor, a continuous decrease of the average duration, due to the experience occurring, so these input values should be dynamically updated.

When a surgery is finished, the name of the surgeon and the time spent are introduced in the DSS. The DSS registers this information in a database and updates the average duration.

*Overview of the DSS.* All the previously explained features are integrated into the DSS whose flowchart is given in Figure 1.

**Figure 1.** Flowchart of the DSS. The steps form the moment when a patient arrives until she/he is surgically operated are shown. The discontinuous colored boxes reference the five main panels included in the developed software tool: (1) Medical teams; (2) Patients; (3) Surgeries; (4) OR timetable; and (5) Scheduling.

When a new patient arrives at the service, she/he is introduced into the waiting list of a medical doctor. Each surgeon is responsible for including their patients in the system. The DSS recognizes the surgeon (using a personal password) and she/he must enter some personal information of the patient, the pathology and the priority of the surgery. Additionally, the DSS saves the actual date (to compute the waiting time in the list) and the surgeon to whom the patient belongs. When a team scheduler wants to perform a scheduling, she/he is recognized by the DSS and automatically composes the ordered waiting list of his patients.

The DSS has a database, which is updated every time that a surgery is performed. The pathology and its surgery durations are saved in the database. Thus, considering this information, the tool assigns average theoretical duration and standard deviation to each surgery in the waiting list. In this way, the vectors *μw* and *σw* are generated. On the other hand, the DSS saves the OR time-table introduced previously by the head of the department—that is, the time blocks booked for each team as well as their duration defined as *l*(*bi*). The DSS performs an operation scheduling in an iterative way. The input data that the team's manager has to introduce in the DSS to schedule the next time blocks are: (i) the minimum confidence level *Cl* and (ii) the number of time blocks to schedule |B|.

When a patient has been scheduled, his/her state changes from "pending" to "scheduled". Once a specific surgery is performed, the surgeon indicates this situation in the DSS and the corresponding patient is removed from the waiting list. Additionally, the surgeon introduces the operating time in the tool and these new input data are used to update the average duration. If finally a scheduled surgery is not performed, the DSS changes the state of this surgery from "scheduled" to "pending".

Based on the previously described DSS, a Software tool for scheduling elective patients in Spanish Hospitals called CIPLAN has been developed. A first version of this software was described in [12]. CIPLAN has been developed in JAVA language and has a SQL database. Moreover, in order to provide a safe work environment and a better security of the medical data, different access levels are provided to the users. The preferences and concerns of medical doctors in the OSD have been considered not only at the surgery scheduling level but also at the user interface level. Thus, a friendly interface divided into six main panels was proposed. These panels are based on the DSS and the first five are shown with color boxes in the flowchart of Figure 1: (1) *Medical teams*, (2) *Patients*, (3) *Surgeries*, (4) *OR timetable*, (5) *Scheduling*, and (6) *User*. Figure 2 shows the interface of the scheduling panel.

**Figure 2.** Scheduling panel. The main interface of the software tool is shown. Particularly, the scheduling panel where once the number of blocks to schedule, the minimum confidence level, and the surgical team are introduced, a summary of the scheduling performed is shown at the bottom part.

#### **5. Simulations**

In this section, some scheduling results derived of simulation are described and analyzed. The approach presented here is compared with two approaches in literature.

#### *5.1. Simulation Methodology*

In order to test the SHA and to compare it with the approaches proposed in [19,20], a discrete event simulation model of the scheduling has been implemented. It is used to simulate scheduling decision for a medical team in the Orthopedic Surgery Department (OSD) at the HCU. One year and the half length (78 weeks) is set for each simulation run. The new patients are assumed to arrive according to a Poisson distribution with an average of six per week. We assume initially a waiting list of 100 patients and every week the new arrived patients are added to the list.

During the last two years, considering historical data in the OSD, the occurrence probability of each surgery type and their duration have been computed. These values are considered in the simulation for generating the surgery of the arrival patients. Figure 3 shows the occurrence of the most common surgeries performed in the OSD while Table 3 indicates their average and standard deviation. Notice that seven types of surgeries represent around 73% of the total surgeries performed; however, all of them are used in the simulations.

**Figure 3.** Occurrence of surgeries. The occurrence of some types of surgery that can be performed in the OSD of the LBH are shown.

**Table 3.** Duration of the surgeries performed in the OSD of the LBH.

For each set of simulations, 50 replications are performed. Each replication is a schedule of 156 time blocks (78 weeks × 2 time blocks per week) all of them having the same time duration of *l*(*bi*) = 390 min (6.5 h). The average values of the probability of not exceeding the total time (confidence level) and the occupation rate are computed. Moreover, the total overtime and the total number of treated patients per year are considered.

In order to be able to compare two different schedules from the point of view of the order of the patients, indicator Ω is defined. This indicator measures the disorder of the patients in the obtained scheduling, so, the smaller it is, the more ordered are the patients. To compute this value, for each time block, an interval [ *fi*, *li*] is defined. If the preference order of the surgeries scheduled in the time block *bi* belongs to the interval [ *fi*, *li*], Ω is not increased. On the contrary, each patient *wj* with a preference order *j* outside the interval increases the value of Ω. The formal calculation of Ω is given in Algorithm 1, where *N p* is the total number of patients scheduled and *Pd* = *N p* |B| is the average number of patients scheduled per time block.

Algorithm 1 computes first the lower (*fi*) and upper (*li*) preference orders of the patients that should be ideally scheduled in the time block *bi*. However, we allow one block before and one block after without any penalization and, for this reason, *fi* = *Pd* · (*i* − 2) and *li* = *Pd* · (*i* + 1). Then, for all patients scheduled in *bi*, if their preference order in the waiting list does not belong to [ *fi*, *li*], Ω parameter is increased with a value depending on the deviation (step 4).

**Algorithm 1** Calculation of Ω parameter in a scheduling of |B| time blocks

**Input:** B, *Pd* **Output:** Ω **1** Ω := 0 **forall** *bi* ∈ B **do 2** *fi* := (1, *Pd* · (*i* − 2)); *li* := *Pd* · (*i* + 1); **forall** *patient (wj) scheduled in block bi* **do 3 if** *the preference order (j)* ∈/ [ *fi*, *li*] **then 4** Ω := Ω+min(|*j* − *fi*|, |*j* − *li*|)

#### *5.2. Simulation Results*

Table 4 shows a comparison between the scheduling obtained by using: (1) the first-fit probabilistic rule (commonly used in hospitals), (2) the batch scheduling approach in [19], (3) the constructive algorithm proposed in [20], and (4) the SHA approach presented in this paper. The scheduling is obtained following the simulation methodology explained in the previous subsection for a minimum confidence level *Cl* = 70%. The average value of probability of not exceeding the total time, occupation rate, and the order of the patients (Ω) are analyzed. Moreover, the total overtime and the total number of treated patients per year are shown.

**Table 4.** Comparison of the one and the half years scheduling using different chain-constrained approaches with a minimum confidence level of 70%.

The three approaches analyzed in Table 4 improve the occupation rates obtained by using the first-fit probabilistic rule. However, the improvement in the occupation rate of the time block implies a decrement in the confidence level. For example, the Batch Scheduling approach [19] achieves the highest occupation rate (79.06%) and the highest number of treated patients (447.78), and, consequently, the lower confidence level (75%) and the highest total overtime (1183 [min]) are obtained. Taking into account the pairs of values occupation rate and confidence level, the Batch Scheduling approach obtains the better solution with: (1) the highest occupation rate and (2) a confidence level within the allowed.

The SHA approach presented in this paper obtains a little worse occupation rate (78.28%) than the Batch Scheduling (79.06%). However, considering the order of the patients given by the indicator Ω, it can be checked that the SHA approach by far obtains the best scheduling (Ω = 395). Medical doctors of the OSD consider that the scheduling obtained using the SHA approach is the best one since it maintains the order of the waiting list better. Due to patients having more orders with SHA scheduling, it is possible to provide an estimation period for the patients in the list and this is a possible extension.

#### **6. Case Study Results**

This section presents a comparison between a scheduling obtained manually in the OSD and the scheduling obtained by using CIPLAN.

#### *6.1. Case Study Description*

**Hospital Description.** The HCU is a public hospital located in Zaragoza, Spain providing health services to around 325,000 people, and it is a reference center of specialized attention of a population over 1 million inhabitants. The hospital has 800 beds and 15 ORs for major surgery, which is performed by eight different surgical departments. The case study considers the Orthopedic Surgical Department (OSD).

The OSD is composed of five medical teams and has assigned 3 ORs (*OR*1, *OR*2 and *OR*3) from 8:30 a.m. to 3:00 p.m. every day. Although there are other time blocks available for the department, the case study considers only these three ORs because the other ones have a variable time-table depending on the urgent surgeries. *OR*1 and *OR*2 are used to perform surgeries of elective patients, while OR3 is assigned to perform urgent surgeries. Table 5 shows an example of the weekly OR time-table of this OSD.

**Table 5.** OR time-table in the Orthopedic Surgical Department of the HCU.

According to Table 5, each medical team has two time blocks per week for performing surgeries on the elective patients.

**Fri** team 3 team 4 team 5

**Methodology of the case study.** The main objective of the developed software tool is to improve the efficiency and quality of the surgical service. Thus, in the case study, the scheduling obtained manually is compared with the scheduling obtained by using CIPLAN. The methodology is as follows:

1. Analysis of manual scheduling. Considering historical data of one medical team of the OSD, 36 consecutive OR time blocks of elective patients are analyzed. For each time block, the following information is obtained:

The effective time in a surgery block is the time of using the OR for surgeries; therefore, the occupation rate in a time block is computed as the effective time divided by the total time.

2. Create the waiting list. The manual scheduling and consequently the surgery types in the waiting list are known. However, the preference order of each surgery in the waiting list is unknown. In order to be able to obtain an ordered waiting list from which the scheduling using CIPLAN is performed, it is assumed that the position of the patients (surgeries) in the list corresponds to the order in which they were operated by using the manual scheduling. That is, the first patient in the waiting list is the first patient operated in the first time block. Thus, considering the surgeries performed in the 36 time blocks, an ordered waiting list containing 111 patients is created.

3. Scheduling using CIPLAN. The data necessary to perform the scheduling using CIPLAN are included in the tool and then the scheduling is performed:

4. Analysis of the scheduling performed using CIPLAN. For each time block scheduled, the following parameters are obtained:

The scheduling obtained by using CIPLAN has not been applied in practice. Thus, to compute the real occupation rate, the time spent when the surgeries were performed (manual scheduling) has been considered. Moreover, in order to obtain the real ending time, 20 min of cleaning time between two consecutive surgeries and 10 min of delay with respect the starting time have been considered.

#### *6.2. Real Case Study*

Figure 4a shows the *expected occupation rate* in each time block. The average value obtained using CIPLAN (77.44%) is greater than using the manual method (75%). Moreover, in the scheduling obtained using CIPLAN, the values of the expected occupation rate in each time block are more concentrated around the average value (*σ* = 3.61) than using the manual scheduling (*σ* = 7.42).

**Figure 4.** Comparison of manual and CIPLAN scheduling. A set of 24 time blocks has been scheduled manually and using CIPLAN. The following values are shown: (**a**) expected occupation rate; (**b**) real occupation rate; (**c**) confidence level of not exceeded the available time; (**d**) ending time.

In Figure 4b, the *real occupation rate* obtained in each time block is shown. Again, the average value obtained using CIPLAN (77.86%) is greater than the average value using the manual method (75.03%). Using CIPLAN, the improvement in the occupation rate with respect to use the manual method is 2.83%.

In Figure 4b, it is also possible to see that using CIPLAN there are some time blocks (1, 2, 23, and 24) with occupation rates too high (about 90%) and could lead to exceeding the available time. This does not happen in manual scheduling. Analyzing the results, it has been observed that it happens because the scheduler (in the manual method) knows the surgery skill of the surgeon and can estimate better the expected duration of the surgery depending on the surgeon. However, the CIPLAN scheduling algorithm initially uses the average duration of the surgeries. These average durations are computed considering historical data, independently from the surgeon. As the tool is used, the time spent in each surgery as well as the surgeon who performs it will be registered in CIPLAN. Thus, the surgery duration will be more and more customized depending on the surgeon and better results will be obtained.

In Figure 4c, the confidence level of not exceeding the available time is shown. The corresponding values are computed by considering the average duration and the standard deviation of each surgery. It can be seen that using CIPLAN does not exist time blocks with a confidence level lower than 69%. However, using the manual method, there are several time blocks with a really low confidence level. For example, in time block 4, the confidence level of not exceeding the available time is only 23.3%. Again, these low values are obtained since the scheduler knows the surgery skills of the surgeon. Particularly, in time block 4, three knee arthroplasties have been scheduled. Considering the average duration and the standard deviation based on historical data, the expected occupation rate is 94.8% and the confidence level that will not exceed the available time is 23.3%. However, the scheduler knows that the surgeon who has to perform the knee arthroplasties needs 25 min less than the average duration of the team. Consequently, the real occupation rate obtained is 73.3% and the time block ends at 13:52.

Figure 4d shows the ending time of each block. It can be observed that, using CIPLAN, there are four time blocks (1, 2, 13, and 24) ending around 3:00 p.m. and the time block 23 exceeding the available time 30 min (ending 15:30). This situation does not happen in manual scheduling.

According to the results obtained in the case study, the quality of the service decreases, because some blocks exceed the available time. However, this situation will be solved when the average duration and standard deviation of the surgeries will be customized depending on the surgeon. This customization will allow that the real occupation rate will be more similar to the expected one and, consequently, not only the number of blocks exceeding the available time will decrease, but also the total exceeding time.

In Figure 5, the time execution of the manual (a) and automatic (b) scheduling are shown by two bar-graphs.

Each row is composed by boxes that represent different actions in one surgical time block. Colored boxes represent surgeries while white boxes are idle times. The idle time could be: (i) delay regarding the starting time if the box is at the beginning or (ii) cleaning time if the box is between two surgeries. On the other hand, to check how much the order of patients in the waiting list is respected, the preference order of the patients inside the colored boxes has been included. Since the waiting list of patients has been constructed considering real surgeries, in Figure 5a, the patients are perfectly ordered. However, considering the automatic method (Figure 5b), although the SHA tries to make the scheduling respecting as much as possible, the order of the patients in the waiting list results in being a little disorderly. In order to see it visually, three colors have been used with a different meaning:

**Figure 5.** Comparison of the time execution of the surgical blocks. The time execution of the 24 time blocks scheduled (**a**) manually and (**b**) using CIPLAN are compared. Numbers inside the boxes represent the preference order of the patients. In case (**b**) color green/orange/red means that the patient has been advanced or delayed at most 1/2/3 time blocks with respect to the manual scheduling.

Remember that the waiting list has been constructed considering the real surgeries performed in 36 surgical blocks, obtaining a waiting list composed of 111 patients. Considering the manual scheduling, patients with preference order from 1 to 73 have been scheduled in the first 24 time blocks. Using CIPLAN, patients with preference order from 1 to 76 (except 65 and 71) have been scheduled in 24 time blocks. The number of patients advanced or delayed two time blocks (orange) is 14 while the number of patients advanced or delayed three time blocks (red) is 9. However only three patients (15, 18, and 57) have been delayed three time blocks and seven patients (22, 25, 47, 54, 58, 60, and 64) have been delayed two time blocks. This means:

#### **7. Discussion**

#### *7.1. The Use and Experiences with the Scheduling Software Tool*

Medical doctors who have used CIPLAN consider it an intuitive, rapid, and efficient software solution that can make easier the scheduling task. With the manual method, the confidence level of not exceeding the available time is unknown before building the schedule. However, with CIPLAN, they find an important advantage that the input parameter is the confidence level and the tool maximizes the occupation rate respecting this confidence level. Moreover, they said that the time spent for scheduling is drastically reduced, a check of the obtained solution being necessary.

Using CIPLAN, from an economic point of view, the improvement in the occupation rate using the manual method is 2.83%. Only considering the OSD, this 2.83% implies an increment Δ*Oc* of 95.65 effective hours per year in the use of the ORs:

$$\begin{aligned} \Delta\_{Oc} &= 0.0283 \cdot 6.5 \left[ \frac{hours}{block} \right] \cdot 2 \left[ \frac{block}{team} \right] \\ &\cdot 52 \left[ \frac{week}{year} \right] \cdot 5 \left[ team \right] = 95.65 \left[ \frac{efe.hours}{year} \right] \end{aligned} \tag{8}$$

Assuming an occupation rate of 77.86% and time blocks of 6.5 h, these 95.65 effective hours are equivalent to an increment Δ*Tb* of 18.89 blocks per year:

$$\Delta\_{Tb} = \frac{95.65 \left[ \frac{eff \,\mu\text{-}h}{year} \right]}{0.7786 \left[ \frac{eff \,\text{-}h}{\text{h}} \right] \cdot 6.5 \left[ \frac{h}{Hoch} \right]} = 18.89 \left[ \frac{blocks}{year} \right] \tag{9}$$

Considering that the cost to open an OR during 6.5 h is 5850 euros, the use of CIPLAN in the surgical scheduling could suppose a saving of about 110,000 euros per year. The savings could be really greater extending the use of CIPLAN to other Departments.

#### *7.2. Future Work*

Doctors propose two new features that could improve significantly the tool. The first one is a new panel to give the possibility of performing virtual scheduling in which all patients included in a waiting list are scheduled. These provisional scheduling will allow doctors to give patients an estimated date of surgery from the moment when he/she is included in the waiting list. On the other hand, the virtual scheduling allows the head of the department to know the number of time blocks necessary to schedule all patients in the waiting list. Thus, he/she will be able to demand from the health administration the exact blocks needed to operate all patients in a certain period of time.

The second feature proposed is to include a new panel in the tool that, once a time block has finished and the starting and ending time of each surgery performed has been introduced in the tool; a comparison between the expected scheduling and the real execution of the time block is shown. In this way, it will be possible to detect the causes of time deviation between the expected and the real durations.

#### **8. Conclusions**

In Spanish hospitals, doctors are responsible for scheduling of elective patients. Currently, this task is performed manually and has three main problems: (i) doctors need to spend their time in the administrative task (scheduling); (ii) usually under or overutilization of the ORs is obtained; and (iii) objectivity can be questioned due to fact that the scheduling is done by a human. This paper presents a *Decision Support System* (DSS) that uses a *specific heuristic algorithm* SHA, in order to help the hospital managers in the scheduling of the elective patients.

Considering the duration of the surgeries and of the cleaning time as random variables with normal probability density function, the SHA with some statistics concepts can be used to schedule the elective patients from the waiting list to the next time blocks in such a way that: (i) the occupation rate of the time block is maximized; (ii) ensures a minimum confidence level of not exceeding the available time; and (iii) respect as much as possible the order of the patients in the waiting list. In addition, the DSS also includes features enabling the management of the medical teams and the OR time-table, updating the waiting list of patients, iterative planning, and automatic improvement of the input data.

In collaboration with the Orthopedic Surgery Department of the "*Lozano Blesa*" Hospital in Zaragoza, a software tool (CIPLAN) based on the proposed DSS has been developed. The preferences and concerns of doctors have been considered not only at the surgery scheduling level, but also at the user interface level. Thus, a friendly interface divided into six main panels is proposed in CIPLAN. It has been developed in Java language and uses a SQL database. Moreover, in order to provide a safe environment to work and better security of the medical data, different access levels are provided to every user.

Finally, in order to check the impact of the software tool in the efficiency and the quality of the surgical service, a case study considering real data from the Orthopedic Department of the HCU has been discussed. The results show that, using CIPLAN, an improvement in the occupation rate of the ORs can be reached.

However, it has been observed that there are significant differences in the surgery duration depending on the experience and on the surgery skill of the surgeon. Thus, if the SHA works with general average duration, over and under occupation rates are obtained more frequently than if it works with customized data. It would be convenient to customize the data of the surgical durations depending on each specific surgeon. In this way, the expected occupation rate and the expected ending time will be more similar to the real ones improving the quality of the service.

As future work, a new panel in CIPLAN can be developed that, once a time block has been finished, allows doctors to compare in a visual way the expected scheduling and the real execution. On the other hand, a new feature in CIPLAN can be included that enables them to perform virtual scheduling in which all patients in the waiting list are scheduled. In this way, an approximated date could be given to patients from the moment when they are included in the system. Moreover, other extensions can be developed—among others, an app connected to CIPLAN allowing patients to know the approximated date of their surgeries.

#### **9. Patents**

The developed software tool has been registered by the Office of Transfer of Research Result (OTRI) at the University of Zaragoza.

**Author Contributions:** All authors contributed to the design of the mathematical model and the conceptual design of the decision support system. J.A. was responsible for collecting the necessary surgical data to analyze the mathematical models and perform the case study. D.C., C.M. and M.S. were involved in the development of the proposed heuristic method as well as the implementation of the DSS. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work has been partially supported by CICYT-FEDER (Spain-EU) under Grant DPI2014-57252-R and by the Aragonese Government (Spain).

**Acknowledgments:** We would like to acknowledge with much appreciation Lidia Castan's help in the collection of the necessary surgical data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **A Prostate MRI Segmentation Tool Based on Active Contour Models Using a Gradient Vector Flow**

**Joaquín Rodríguez 1, Gilberto Ochoa-Ruiz 2 and Christian Mata 3,\***

Received: 28 July 2020; Accepted: 2 September 2020; Published: 4 September 2020

**Abstract:** Medical support systems used to assist in the diagnosis of prostate lesions generally related to prostate segmentation is one of the majors focus of interest in recent literature. The main problem encountered in the diagnosis of a prostate study is the localization of a Regions of Interest (ROI) containing a tumor tissue. In this paper, a new GUI tool based on a semi-automatic prostate segmentation is presented. The main rationale behind this tool and the focus of this article is facilitate the time consuming segmentation process used for annotating images in the clinical practice, enabling the radiologists to use novel and easy to use semi-automatic segmentation techniques instead of manual segmentation. In this work, a detailed specification of the proposed segmentation algorithm using an Active Contour Models (ACM) aided with a Gradient Vector Flow (GVF) component is defined. The purpose is to help the manual segmentation process of the main ROIs of prostate gland zones. Finally, an experimental case of use and a discussion part of the results are presented.

**Keywords:** Active Contour Models; snake segmentation; GVF; prostate imaging

#### **1. Introduction**

Prostate cancer (PCa) remains one of the most commonly diagnosed solid organ tumor types among men. Early diagnosis and an active follow-up can allow improved prognosis and prevent life-threatening conditions. Once the decision of treatment has been taken, having the most complete information for treatment and follow up is crucial. In 2020, there will be approximately 1,806,590 cancer cases diagnosed and 606,520 cancer deaths are projected to occur in the United States. This classification has been reported by the ICD-O, excepting childhood and adolescent cancers, which were classified according to the ICCC [1]. On the other hand, PCa is the third predicted cause of death in EU with 78,800 death in 2020 [2]. This entails a rate estimation of 10.0/100,000 men.

Medical support systems used to assist in the diagnosis of prostate lesions generally focus on prostate segmentation [3–6]. They rely on computerized techniques for prostate cancer detection applied to ultrasound, magnetic resonance and computed tomodensitometric images [7]. For example, Vos et al. [8] used 3D-T2 imaging to define a specific Region of Interest (ROI), which was subsequently used on diffusion- and perfusion-weighted images to extract relevant features. MRI has been established as the best imaging modality for the detection, localization and staging of PCa on account of its high resolution and excellent spontaneous contrast of soft tissues and the possibilities of multi-planar and multi-parametric scanning [9]. Among the techniques used to detect PCa, MRI also allows non-invasive analysis of the anatomy and the metabolism in the entire prostate gland.

Furthermore, many recent works confirm that the combination of several MRI techniques facilitates the evaluation of prostate diagnoses according to the PI-RADS classification, both for radiologists and medical experts [10]. For instance, a thorough evaluation of various prostate cancer studies was presented in [11]; this study highlights that the correct localization of ROI-containing tumor tissue as the most prevalent problem encountered in the diagnosis of a prostate cancer and other related tasks. Traditionally, clinical experts use different software tools to validate the diagnoses and make annotations in several files. This approach is far from being efficient when managing and dealing with abundant medical data. For the reasons stated above, we propose a new software tool to provide support for the time consuming segmentation process used for annotating images in the clinical practice. Additionally, it enables radiologists to make use of novel semi-automatic segmentation techniques instead of the traditional manual segmentation approaches.

The contributions of this paper are two-fold: we introduce a novel prostate cancer segmentation graphical user interface on which several segmentation algorithms can be easily implemented. In order to demonstrate the strengths and the potential applications of the proposed tool, we included several segmentation algorithms with the software:

The extensions of original algorithm, in tandem with the novel graphical user interface greatly simplify the segmentation process as it enables a boosted semi-automatic interaction with the clinicians who, eventually, can validate the correctness of the ROI generated by the different methods. Finally, a quantitative study to validate the correctness of the generated segmentation compared to the ground truth and the algorithms is provided. It is worth to mention that the developed program serves as an initial version of a tool that we hope, will create a lot of opportunities to contribute with newer algorithms or improvements to the existing ones in the field of prostate segmentation.

The rest of this paper is organized as follows: in Section 2 we introduce the used dataset and we explain the implementation and features of the proposed Graphical User Interface. Then we describe the Basic ACS and the modifications we introduced to make it work in prostate images. Finally, we present details about the Gradient Vector Flow (GVF) model and its implementation into our tool. In Section 3, we include a single experimental case study to demonstrate the advantages of the proposed approach. In Section 4 we discuss the results in the light of the clinical perspective, and in Section 5 we present the conclusions as well as future works.

#### **2. Materials and Methods**

In this section we present the source database and the main contributions of the article, namely the proposed GUI for prostate segmentation process and later on, the integration of the Basic ACS model and its extensions using GVF, the GVF ACS model.

#### *2.1. Database and ROIs References*

A database with prostate MRI, based on clinical data with tumor and healthy cases was used. The examinations used in our study contained three-dimensional T2-weighted fast spin-echo (TR/TE/ETL: 3600 ms/143 ms/109, slice thickness: 1.25 mm) images acquired with sub-millimetric pixel resolution in an oblique axial plane. Each study comprises a set of 64 slice images. This database was used in previous works [11] and it has been facilitated for testing the software.

The ground truth annotations were performed using the ProstateAnalyzer tool [11], that allows the drawing of an annotation in a given MRI. For this task, experts with more than 10 years of experience on prostate imaging were asked to provide the annotations using the tool. All the dataset and ground truth data are provided from the medical imaging department of the University Hospital of Dijon (France). The institutional committee on human research approved the study, with a waiver for the requirement for written consent, because MRI and MRSI were included in the workup procedure for all patients referred for brachytherapy or radiotherapy.

The prostate is composed of a peripheral zone (PZ), central zone (CZ), transitional zone (TZ) and an anterior fibromuscular tissue (AFT) depicted in Figure 1. Most cancer lesions occur in the peripheral zone of the gland, fewer occur in the TZ and almost none arise in the CZ. A detailed description of the influence of the prevalence factor risk according the prostate zone is described in [13]. The T2WI modality was chosen because it provides the best depiction of the prostate's zonal anatomy.

**Figure 1.** Anatomy of the prostate in (**a**) transversal and (**b**) sagittal planes [11].

The PZ represents up to 70% of a normal prostate gland and around 75% of prostate cancers originate in this zone. The CZ represents about 25% of a normal healthy prostate gland in young adults. Even though the frequency of cancers originating here is much lower, they tend to be of the more aggressive type [14]. An example of a real study case provided by the panel visualization of the prostate tool [11] is depicted in Figure 2. All the overlays are annotated with manual drawings made by radiologists. Each region is represented by different colors (CZ in white and PZ in blue) focusing on a tumor area represented in red (Tum). This example shows the importance of making annotations in regions of interest. In this case, a manual segmentation is validated by experts.

**Figure 2.** An example of the anatomy (T2WI) panel visualization of the tool developed in [11] with the ROIs annotated by radiologists.

Therefore, one of the main motivations of this work is to develop a new tool focused on improving the results of this initial manual segmentation using Active Contour Models (ACM), facilitating the task of the experts.

In order to attain the above-mentioned goals, we have implemented a novel Graphical User Interface that permits to import, analyze and annotate prostate MR images. The program requires an initial annotation from the user, that does not require to be perfect, and the outcome is an annotation more accurate than the initial one. The modular design of the software allows for integrating segmentation methods as plug-ins or extensions that can be used for testing and clinical purposes. Additionally, this design enables easy performance comparisons between already integrated methods and new methods. In what follows, we will discuss the design of this GUI and subsequently, we will present a study case using three implementations of an ACS model for segmenting the different prostate sections in MR images. It is important to note that this tool and the developed algorithms have not yet been clinically tested and the results presented here represent a proof of concept. However, for the design of the GUI, the expertise and suggestions of clinical experts have been taken into account.

#### *2.2. Implementation of the Graphical User Interface for Prostate Segmentation*

In this section, the code organization and the user interface of the program will be explained. The main window of the program is shown in Figure 3. All the code has been developed in Python, using the modules PyQt for the user interface and PyDicom to open the files included in the testing dataset. The code has been organized in modules that do not interact with each other. The integration of all of them has been made in the Main Window class.

**Figure 3.** GUI interface of the developed active contour segmentation tool using a prostate MR image.

Each module is placed in folders as follows:

• **DICOM.** This module contains the code that is in charge of handling the Digital Imaging and Communication On Medicine (DICOM) formatted files. It can open either a folder or a single DICOM image. Each file is assumed to contain only one image. The read information will be stored locally on RAM, and after the file or files has/have been opened, a signal will be emitted, and the information can be requested as many times as the user wants.

The data available to the user includes the images and the patient data fields. This module is also in charge of anonymizing the files and writing them into disc. The functionalities of this module can be accessed through the buttons in the Section 1 of the Figure 3.

• **Image segmentation.** This module stores all the Python code related to the segmentation procedure. The main class in this module is Segmentation-Hub, which is in charge of selecting the right class for the object that will produce the segmentation, according to the method chosen by the user. The selection of the segmentation method can be done in Section 2 of Figure 3. By implementing the methods in this manner, the actual implementation of the segmentation task is completely transparent for other classes that make use of this module. To select certain segmentation method, only a label that defines the method must be provided.

This module also contains two other files with useful functions:

Each button will request the main window to ask the user to enter the initial contour in the left image. Once the user has finished, the gathered contours are returned back to this module.

Afterwards, the corresponding masks for those contours can be generated for the snake algorithm (the mask consists of an image with values of 1 in the interior of the contour, and zeros outside of it). The interface for this module is represented in Section 6 of Figure 3.

• **User data widget.** It is implemented as an efficient and easy to use manner to handle all the user information as an object, and to show it. This class stores all the information related to the patient or the study, retrieved from the DICOM file. It also contains the corresponding widgets that are shown below the left image area, where the gathered information is displayed (see Section 7 of Figure 3).

A block diagram of the proposed tool components is depicted on Figure 4.

**Figure 4.** Main software components of the proposed tool. The Segmentation-Hub module can be extended with new modules seamlessly.

Finally, regarding the user interface showed in Figure 3, this program has four well-defined sections:

To the best of our knowledge, tools as the proposed in this work, which integrate a Graphical User Interface and semi-automated segmentation algorithms for image annotation do not exist (at least not for prostate MR images). Therefore, a preliminary effort to develop a GUI that includes several active contour model implementations is very important and has driven the development of this work. On the other hand, from the implementation point of view, this program represents an excellent tool to deploy new snake-based algorithms: given that the initial contour and the original image are available in the Segmentation-Hub module, new algorithms can be created using Python, and added as an standalone file to the Image Segmentation folder. This procedure will add another option to the GUI, and the new algorithm's performance and accuracy can be compared with the already existing techniques. This code-base is now available as a GitHub repository to the community [15].

In order to demonstrate how different segmentation methods can be integrated into the proposed tool, we have implemented various versions of an Active Counter Model within the Image Segmentation module of our tool. These methods have been included not only for demonstration purposes, but also to serve as a baseline for future developments in prostate MRI segmentation research (people can improve them, or develop new methods based on them). The proposed Active Contour Model, to be discussed in the following section, also represents a novel contribution, as it has not been previously reported to be used for prostate segmentation purposes.

#### *2.3. Implementation of the Proposed Active Contour Model*

The base method implemented is inspired on [12,16]. Both of them are active contours—level-set algorithms. Particularly, the model by Pluempitiwiriyawej et al. [12] was conceived and customized to work with MR images. The general idea of snake algorithms is to modify a curve based on certain forces that push the boundaries of such curve until it reaches the border of the object of interest. A particular feature of this type of segmentation algorithm is that they are unsupervised methods, given that there is not feedback on whether the boundary is following the correct borders or not. Due to this assumption of snake-based methods, an initial contour must be provided, and the closer this initial contour is to the real object, the better the result of the segmentation.

Typically, most snake-based algorithms follow the information about the gradient of the image, which means the curve will try to evolve towards the edges of the object. Nonetheless, if the object does not contain well defined edges, or it has been previously smoothed in order to reduce the noise, this algorithm might fail in following the shape of the object. As a consequence, the outcome will be either a contour that is smaller or larger than the object's boundaries, or in the worst case, the contour will be reduced until it disappears. This is the reason why several researchers developed specific algorithms depending on the type of object to detect, and the type of image to deal with. This is done adding different terms called *priors*. Each constraint added to the algorithm defines a prior, and its corresponding weight will define its contribution to the force.

In this previous work [12], the goal is to segment the shape of the left ventricle of the heart, but since it is aimed to work with MR images, and the priors are similar to the assumptions that can be done for the prostate zones, we considered it as a good starting point for carrying out our experiments and implement the GVF extension to demonstrate its viability in this image modality. The basic assumption of the developed algorithm is that the contour will split the image into two areas, the target object and the background. The evolution of the curve will be driven trying to minimize an energy functional that in this case, it is composed of four components. Each of them will define the priors the final contour must comply with [12]:

Then the functional that defines the energy described above, depending on the contour *C*, can be expressed as:

$$J(\mathbb{C}) = \lambda\_1 I\_1(\mathbb{C}) + \lambda\_2 I\_2(\mathbb{C}) + \lambda\_3 I\_3(\mathbb{C}) + \lambda\_4 I\_4(\mathbb{C}) \tag{1}$$

where *λi* are the weights of each prior included in the energy function. This function will be used to make evolve the contour. It follows that it is possible to transform the information of the contour in a level set function, i.e., the information about the contours can be expressed as a function, where the pixels that belong to the object have certain value (or level), the background ones have another level, and the frontier, which is the contour *C*, is represented by another level. In the case of this paper, the contour is expressed through the level set function as in Equation (2). As it can be seen, the contour is in the level where the function *φ* is equal to zero. Ω represents the space of the image.

$$\mathbb{C} = \{(\mathbf{x}, \mathbf{y}) \in \Omega : \phi(\mathbf{x}, \mathbf{y}) = 0\} \tag{2}$$

Then, the energy functional can be expressed as depending on *φ*.

$$J(\phi) = \lambda\_1 I\_1(\phi) + \lambda\_2 I\_2(\phi) + \lambda\_3 I\_3(\phi) + \lambda\_4 I\_4(\phi) \tag{3}$$

This function *φ* evolves with each iteration (i.e.,time), depending on the direction in which the forces push the contour. In the following sections the particularities of each term will be described.

#### 2.3.1. Region-Based Model Matching—Term *J*1(*φ*)

The aim of this term is to split the image space Ω into two regions that contain homogeneous pixels: the object and the background. Let's call them Ω1 and Ω2, respectively. These two regions are separated by the contour *C*. An homogeneous region implies a region where all the pixels that are contained in it are described by the same stochastic model. The points that belong to the object have intensities *u*1 = *ui*,*j* : (*i*, *j*) ∈ Ω1 , and they are described by the model *M*1, and the ones that belong to the background have values *u*2 = *ui*,*j* : (*i*, *j*) ∈ Ω2 , and they are represented by the model *M*2.

The idea of the model is as follows: let's consider the image *u*0 to be formed by two regions, of piecewise-constant intensities *c*1 and *c*2, and that the intensities of the object are represented by

*c*1. The boundary of this region is *C*0. Then, the pixel values of the image *u*0 are *c*1 inside *C*0, and *c*2 outside of it. Given this, it is possible to build the following fitting term:

$$F(\mathbb{C}) = F\_1(\mathbb{C}) + F\_2(\mathbb{C}) = \int\_{\text{inside}(\mathbb{C})} |u\_0(\mathbf{x}, y) - c\_1|^2 d\mathbf{x} dy + \int\_{\text{outside}(\mathbb{C})} |u\_0(\mathbf{x}, y) - c\_2|^2 d\mathbf{x} dy \tag{4}$$

where *C* is any variable curve, and *c*1 and *c*2 are two constants that depend on the contour *C*. Based on the previous explanation, if *C* = *C*0, then *F*(*C*) is approximately zero, because each distance between the image value and the constant *ci* are approximately zero for each region. If the curve *C* is outside *C*0, then its interior contains parts of the object and parts of the background, while the exterior of *C* contains only elements of the background. Because of this, *F*1(*C*) > 0 and *F*2(*C*) ≈ 0. In the backwards case, i.e., the contour *C* is smaller than *C*0, the outside contains parts of the object and the background, whilst the inside only contains parts of the object. It follows that *F*1(*C*) ≈ 0 and *F*2(*C*) > 0. Finally, this fitting function can be used to find the boundaries of the object, by minimizing its value such as

$$\inf\_{\mathbf{C}} \left\{ F\_1(\mathbf{C}) + F\_2(\mathbf{C}) \right\} \approx 0 \approx F\_1(\mathbf{C}\_0) + F\_2(\mathbf{C}\_0) \tag{5}$$

In order to achieve this boundary, the evolution of the curve *C* is done by moving pixels from the exterior to the interior of the contour, and vice-versa. The first estimation of the energy function for this prior is in Equation (6):

$$J\_0(\mathbb{C}) = p(\boldsymbol{\mu}|\mathbb{C}, M\_1, M\_2) \tag{6}$$

where *p* is the probability density function of the image intensities *u*, given the contour *C*, and the models *M*1 and *M*2. Over these probabilities, it is assumed that each element, object and background, are statistically independent, and they have pdfs *p*1 and *p*2, respectively. Moreover, the pixel intensities are considered to be independent from each other. With these assumptions, and taking the negative of the logarithm of the pdfs, the Equation (6) becomes Equation (7).

$$J\_1(\mathbb{C}) = \sum\_{(i,j)\in\Omega\_1} -\ln(p\_1(u\_{i\bar{j}})) + \sum\_{(i,j)\in\Omega\_2} -\ln(p\_2(u\_{i\bar{j}})) \tag{7}$$

Assuming a continuous function, the contour embedded in the zero level-set function *φ*, and considering the Heaviside function in order to create a single integral over all the domain of the image, the Equation (7) can be generalized as shown in Equation (8).

$$J\_1(\phi) = \int\_{\Omega} -\ln[p\_1(u(\mathbf{x}, y))]H(\phi(\mathbf{x}, y)) - \ln[p\_2(u(\mathbf{x}, y))](1 - H(\phi(\mathbf{x}, y)))d\mathbf{x} dy \tag{8}$$

It is assumed that both elements in the image follow a Gaussian distribution, with means *m*1 and *m*2, and variances *σ*2 1 and *σ*2 2 , which are unknowns, and for each contour *C*, they will be estimated.

#### 2.3.2. Edge-Based Function—Term *J*2(*φ*)

This term will push the contour towards the edges of the object. The computation of this type of force is done through the edge-map *Υ*(*x*, *y*) generated from the image *u*(*x*, *y*). The most common implementation is the gradient of the smoothed image.

$$\mathcal{Y}(\mathbf{x}, \mathbf{y}) = -|\nabla G\_{\sigma} \ast \mathbf{u}(\mathbf{x}, \mathbf{y})| \tag{9}$$

where *Gσ* is a Gaussian kernel of parameter *σ*, ∇ is the gradient operator and (.) ∗ (.) is the convolution operator. Finally, the energy functional, depending on the level set function *φ*, which corresponds to this prior is given by the Equation (10).

*Appl. Sci.* **2020**, *10*, 6163

$$J\_2(\phi) = \int\_{\Omega} \mathcal{Y}(\mathbf{x}, y) \delta(\phi(\mathbf{x}, y)) |\nabla \phi(\mathbf{x}, y)| dx dy \tag{10}$$

where *δ*(*φ*(*x*, *y*)) will return non-zero values only for the pixels that are on the contour *C*.

#### 2.3.3. Object Shape Prior—Term *J*3(*φ*)

This term has been introduced in [12] to reinforce an elliptical shape in heart MR images, but it can be changed by any other shape model depending on the specific object to be segmented. The addition of this prior will improve the performance of the general active contour implementation. For instance, if the target to segment has a circle-like shape, then the snake function can be constrained to satisfy an ellipse equation. In that case, the functional corresponding to this energy term can be defined as *J*3(*φ*) = *C D*2(*x*, *y*)*ds* where *D*(*x*, *y*) is the ellipse distance function defined as:

$$D(\mathbf{x}, \mathbf{y}) = a\mathbf{x}^2 + b\mathbf{x}\mathbf{y} + c\mathbf{y}^2 + d\mathbf{x} + e\mathbf{y} + f \tag{11}$$

The parameters *θ* = [*a*, *b*, *c*, *d*,*e*, *f* ] are computed in such a way that the cloud of points around the contour *C* fits with this ellipse equation. This means that the parameters of *θ* must be computed for each contour *C*, and they will change the force direction in order to minimize this distance function. Finally, introducing the level set function *φ*, where the contour *C* is embedded, the shape prior is computed as in Equation (12), when integrating in all the space Ω.

$$J\_3(\boldsymbol{\phi}) = \int\_{\Omega} D^2(\mathbf{x}, \boldsymbol{y}) \delta(\boldsymbol{\phi}(\mathbf{x}, \boldsymbol{y})) |\nabla \boldsymbol{\phi}(\mathbf{x}, \boldsymbol{y})| dx dy \tag{12}$$

#### 2.3.4. Contour Smoothing—Term *J*4(*φ*)

Finally, it is desirable that the generated contour is smooth. This can be achieved if the Euclidean arc length of the contour *C* is minimized. It follows that the contribution function for this prior is the one represented in Equation (13).

$$J\_4(\mathbb{C}) = \int\_{\mathbb{C}} ds \tag{13}$$

It is worth mentioning that if only this term is included in the functional, the shape of the contour will become a circle, and then it will start shrinking until it eventually disappears. Nonetheless, when all the terms are put together, this situation is avoided, and the effect is to smooth the curve. Again, adding the level set function where the contour is embedded, the final function for *J*4 is defined in Equation (14).

$$J\_4(\phi) = \int\_{\Omega} \delta(\phi(\mathbf{x}, y)) |\nabla \phi(\mathbf{x}, y)| dx dy \tag{14}$$

The following step is to combine all the functional together, and expand the minimization problem, in such a way the function *φ*(*x*, *y*), which contains the contour *C* as its level zero, evolves towards the edge of the object. The entire procedure is well explained in [12], and it consists of three tasks:

$$\frac{\partial \phi}{\partial t} = \left[\lambda\_1 (M\_2 - M\_1) - \nabla P \frac{\nabla \phi}{||\nabla \phi||} - P\kappa \right] \delta(\phi(\mathbf{x}, y)) \tag{15}$$

where *M*1, *M*2 are the stochastic models of the object and the background, respectively, *κ* = *div* ∇*φ* ||∇*φ*|| is the curvature of the contour, *φ*(*x*, *y*) is our level set function, and *P*(*x*, *y*) = *λ*2*Υ*(*x*, *y*) + *λ*3*D*2(*x*, *y*) + *λ*4 is the potential function.

The different values of *λ*1, *λ*2, *λ*3 and *λ*4 will define the contribution of the corresponding energy term, from *J*1 to *J*4, and the one with the biggest value will enforce its properties over the others. In practice, the parameters setting is done empirically, and there is no rule that defines how to chose them. It depends on the quality and the properties of the input image. One solution is to make the weights to evolve at each iteration to obtain the best results. The general idea of this evolution is that at the beginning, the algorithm should minimize the energy functional trying to cover all the object, pushing the boundaries towards the edges (effect produced by *λ*1 and *λ*2), and when the algorithm is close to the end, the curve must enforce the shape of the contour and its smoothness (*λ*3 and *λ*4).

For the particular case of prostate segmentation, several objects with different shapes have to be considered. Therefore, the shape prior cannot be only an ellipse, as it happens in the case of the left ventricle of the heart. Thus for this implementation, the value of *λ*3 has been set to zero.

Additionally, through several experiments, we concluded that the contribution produced by the smoothness prior is not noticeable. In fact, it produces small variations in the contour, and the result is that the algorithm takes more time to converge. This is the reason why the value of *λ*4 has been set to zero too. Finally, the suggested evolution of the weights is shown in Figure 5.

**Figure 5.** Evolution of *λi* weights for the energy functional.

#### *2.4. Gradient Vector Flow Approach (GVF)*

When dealing with snakes algorithms, there are two types of implementations: *Parametric Active Contours (PACs)* or *Geometric Active Contours (GACs)*. They differ in how they define the curve function, since equations formulated for the minimization problem stays almost the same. PACs are based on a parametric curve of a shape represented in Equation (16):

$$\mathcal{C}(s) = \{ [\mathbf{x}(s), \mathbf{y}(s)], s \in [0, 1] \}\tag{16}$$

with the constraint that it should be closed, i.e., *C*(0) = *C*(1). Geometric Active Contours, on the other hand, define the curve to be deformed as embedded in a level-set function, generally as its level zero. This last case is the one implemented and described in Section 2.3. GACs have advantages over PACs,

since they are computationally more stable and they are able to change their topology as the function is deformed. Nonetheless, the implemented method is an edge-based method, since it uses the gradient as pushing force, and local intensity information to determine if a pixel belongs to the object or not. Due to this, several problems can appear, such as:

Additionally, the fact of using gradients only, brings another effect related to concavities, that is illustrated in Figure 6.

**Figure 6.** Example of Gradient vector field of a U shape object. (**Left**) Original image. (**Center**) Gradient of the object. (**Right**) Zoom of the vector field in the central concavity of the object.

In this case, a simple object with a narrow concavity, and its corresponding gradient vector field are shown. From it, it is easy to see that:

These last two problems can be solved in the previously explained algorithm with an external force that, not only adds vectors where the gradient is zero, but that also adds a downwards component in the long concavities in such a way that the active contour will be pushed to the interior of this type of areas. This external force is called Gradient Vector Flow (GVF) [17]. Adding this component to the energy described in Section 2.3 will help to not only make the system more robust against initialization effects, but it will improve the outcome by attaching the resulting snake really close to the edges of the object.

The GVF is a dense vector field **v**(*x*, *y*)=[*u*(*x*, *y*), *v*(*x*, *y*)] derived from the original image, and it is the result of minimizing the energy term defined in Equation (17) [17].

$$\varepsilon = \iint \mu (u\_x^2 + u\_y^2 + v\_x^2 + v\_y^2) + |\nabla Y|^2 \left| \mathbf{v} - \nabla Y \right|^2 dxdy \tag{17}$$

where *ux* and *vx* are the partial derivatives of *u* and *v* with respect to the variable *x*, respectively, *uy* and *vy* are the partial derivatives of *u* and *v* with respect to the variable *y*, respectively, and *Υ*(*x*, *y*) is the edge map of the original image. The parameter *μ* is a constant that serves as regularization: it adjusts the balance between the first and the second term of the integral, and it should be adjusted based on the noise. The higher the noise, the higher the value of *μ*. Interpreting the energy term presented in (17), it is possible to see that if the gradient of the edge map is too small, the energy term depends only on the partial derivatives of *u* and *v*. On the other hand, if this gradient is large, the integral is dominated by it, and the value of the GVF that minimizes the energy is **v** = *Υ*, i.e., the gradient vector flow is equal to the gradient of the edge map on those cases. It can be demonstrated that the vector field **v** = (*u*, *v*) that minimizes (17) is the result of the decoupled pair of differential Equations (18) and (19).

$$
\mu \nabla^2 \mu - (\mu - \mathcal{Y}\_x)(\mathcal{Y}\_x^2 + \mathcal{Y}\_y^2) = 0 \tag{18}
$$

$$
\mu \nabla^2 \upsilon - (\upsilon - \mathcal{Y}\_\mathcal{Y})(\mathcal{Y}\_\mathcal{X}^2 + \mathcal{Y}\_\mathcal{Y}^2) = 0 \tag{19}
$$

where ∇2 is the Laplacian operator. Similarly as in Equation (17), it is possible to see that in the regions where the gradient of the edge map is small, the solution of the differential equations for the GVF are the solution of Laplace's equations. Furthermore, this solution is interpolated from the boundaries of the object, and that's the reason why the resulting field, even far away from the object, points to the edges of it. The gradient vector flow of the same image as in Figure 6, is shown in Figure 7. Now, the solution to the Equations (18) and (19) provides vectors in the areas where the gradient is zero, and the vectors in the concavity contains downwards components too.

**Figure 7.** Example of normalized Gradient Vector Flow. (**Left**) GVF vector field for the same image used before. (**Right**) Zoom in the concavity.

Once the GVF has been computed, it is fixed for the entire the image, and what remains is the optimization equation that will guide the deformation of the active contour to the edges of the object, based on this new force. Several researchers have been working on adding gradient vector flow to the traditional snake equations. For the implementation made in Section 2.3, only GACs algorithms can be adapted to be added as an external force to the energy functional defined in the Equation (15), as it has been done in [18,19]. A simplified version of the GVF optimization equation from [19] is used, and defined as:

$$
\gamma \frac{\partial \phi}{\partial t} = \left[ a\kappa - \upsilon . \frac{\nabla \phi}{|\nabla \phi|} \right] |\nabla \phi| \tag{20}
$$

The parameter *γ* regularizes the step size, *φ* is the level-set function, *κ* = *div* ∇*φ* |∇*φ*| is the curvature, **v** is the Gradient Vector Flow vector field, and the parameter *α* regulates the balance between the two forces. The right hand factor of (20) is added directly to Equation (15), giving Equation (21):

$$\frac{\partial \phi}{\partial t} = \left[\lambda\_1 (M\_2 - M\_1) - \nabla P \frac{\nabla \phi}{||\nabla \phi||} - P\kappa + \beta \left( \left[ a\kappa - v. \frac{\nabla \phi}{|\nabla \phi|} \right] |\nabla \phi| \right) \right] \delta(\phi(x, y)) \tag{21}$$

This new segmentation algorithm has been included in the list of options to segment the different prostate zones. After having implemented the *Basic ACS* in Python, the extension with a GVF component is trivial: all the components of this new force are already available in the code except for the vector field **v**. Since this vector field is fixed for a given image, a new class has been created that solves the optimization problem for the GVF. This new piece of code is executed before the main loop of the segmentation algorithm. Since there are no rule for setting both *α* and *γ*, we arbitrarily decided to set them to 1. In order to adjust the contribution of this new force, the *β* factor has been introduced. For the following experiments, this factor has been also set to 1. In the coming sections, the graphical user interface result of the integration of all the modules will be showed. Then, the described algorithms will be tested to show their performance.

#### **3. Experimental Case of Use for Demonstration**

In this section a case of use for demonstration is presented. Previously, the proposed GUI has been explained in Section 2.2. A detailed definition of this algorithm and the introduced improvements have been explained in Sections 2.3 and 2.4. The following study case demonstrates how an user can use the proposed tool to manage and segment a real prostate image.

First, the user must load an image in DICOM format. For this application, the snake initialization is required; an example of the initial contours for a study is shown in Figure 8. In the area of the GUI that corresponds to the initialization of the snake, the three first buttons are used for selecting the different zones of the prostate, whereas the fourth one gives the possibility of selecting a tumor.

Located below these buttons, there are two radio buttons, enabling to choose the shape of the initial contour, which can be circular or manual. The first option corresponds to an initial contour with a shape of a circle whose radius can be specified with the spin box located at the end of this area. The possible values are between 1 and 30 pixels. The manual input enables the user to enter a random polygon as initial contour. If selected, the user just has to make a click with the mouse in the vertices of the desired polygon. There are no limits for the amount of points, or the shape of the polygon.

The steps required to initialize the snake for any of the prostate areas are the following:

**Figure 8.** Example of snake initialization for the four zones.

The user can repeat the same procedure for all the zones, or just some of them. In order to differentiate the various selected areas, each initial contour will be drawn with a different color. There are no constraints about the type of contours to use. All the contours can be manual, all of them can be circles or they can be a combination of them.

The following area of the control panel is used for triggering the segmentation. This section is depicted in Figure 9. After an image has been loaded, at least one initial contour has been entered, then the segmentation process will start after pressing the *Run* button. The segmentation method will be triggered as many times as initial contours have been entered, and they will be executed sequentially. During this process, the user can cancel the segmentation at any moment, pressing the *Cancel* button. Once the algorithm had finished, the results will be shown in the right image. Additionally, the color code of the contours is the same for both, left and right images. For instance, the resulting contour of color blue in the right image has as initial contour the blue snake from the left image.

If some results have been produced, then a 3D representation of them can be visualized if the user presses the button *Show* 3*D model*. This representation consists only of the resulting contours, replicated several times along the Z axis. It is an extra option for users to visualize the main structure of the segmentation in a 3D model. In fact, this option should be improved using better visualization techniques, but it is out of scope of this work for the moment. The amount of repetitions has been set to 50. The axis are represented in real scale, in mm, if the DICOM file contains the information about slice thickness and spatial resolution. If they are not provided, a pop up message will appear reporting the situation, and the X and Y axis will be presented in pixels, and the Z axis will have a size of 10 units. An example of this 3D representation of a result is shown in Figure 10.

Regarding to file management tasks provided to this tool, users can open a DICOM formatted file or a folder that contains several DICOM formatted files. When the *Open file* button is pressed, a pop up window will appear, and it will only allow the user to select a file, not folders. If a file is not in DICOM format, then the program will reject it, and it will not do further actions. If the user wants to open a folder, it should click the *Open f older* button, navigate until the folder it wants, and go into that folder. Without touching any file, users should click the *Choose* button and the system will automatically search for files, and open only the DICOM formatted files. It is important to emphasize that the search system is not recursive. That means, if there are folders inside the selected

one, and those folders contain DICOM files, the program will not open them. It will only open the files that are in the current directory.

**Figure 9.** Example of snake results after the snake segmentation.

**Figure 10.** Example of 3D representation of the segmented prostate zones.

As it can be seen in the graphical user interface, there are three types of algorithms: *Basic ACS*, *Simplified ACS* and *GVF ACS*. The first and the second algorithms are based on [12], but the second one does not follow the mathematical approach described in this previous work. On the other hand, *Basic ACS* implements exactly what the paper states, as explained in Section 2.3. Both algorithms have been included so that there is an example in the code of how to add a new algorithm to this code-base. Additionally, the *Simplified ACS* performs a faster than *Basic ACS* but with less accuracy. Nonetheless, in order to produce improvements to the original paper, the *Basic ACS* was required.

Finally, considering that this application is focused on DICOM format and it can be used for radiologists and/or medical experts, it is important to maintain the private information of the patient. For this reason, an anonymization module is integrated in this tool. If the user wants to anonymize the opened file or files, it should click into the *Anonymize* button. If there is a file/files opened, a pop up window will be shown where the user must specify the destination folder of the anonymized files. The name of the output files will be the same as the source files, but the string \_*anonymized* will be appended to it. For instance, if the original file is called *Image*00001, the output file will be called *Image*00001\_*anonymized*. The fields that are anonymized are: *Patient Name*, *Patient ID* and the

*Patient Date of Birth*. Those fields are not cleared, but changed by some default data: *Patient Name* is changed by *NN*, *Patient ID* is transformed to *A00000000000*, and the date of birth is changed to 01/01/1970.

We believe that the proposed tool can be beneficial for the research community as it enables to integrate various image management tools thoroughly discussed in previous sections, as well as the capability to segment prostate cancer lesions and other areas as mentioned in the introduction. Even though the application includes the option of segmenting four zones only (CZ, TZ, PZ and Tumor), nothing stops the users of trying to segment other areas of their interest (additional buttons can be appended into the GUI in order to select more areas, or the already existing buttons can be used to segment those areas). These features, as well as the capability to integrate new segmentation models seamlessly make it very appealing for research and teaching purposes in clinical settings. In the former case, it can help in speeding up the creation of datasets for machine learning applications. However, in our experience, doctors are reluctant to use fully automated annotation tools and prefer to use interactive semi-automated tools that enable them to verify and contrast the results with other specialists. In order to test the performance of the proposed tool, we performed one experiment to assess the quality of the segmentation algorithm implemented in this work, which will be detailed in the following section.

#### **4. Discussion of the Segmentation Method Results**

In this section, we present a deep analysis of the results of a single segmentation experiment performed using the implemented algorithms, namely the *Basic ACS*, *Simplified ACS* and *GVF ACS* versions of the snake-based models. For all the cases, the initial contours used for all the zones are shown in Figure 11.

**Figure 11.** Initialization used for the testing of all the algorithms.

The results for all the algorithms and the ground truth for the selected image are shown in the Figure 12.

From this experiment, it is possible to confirm what it has been said in previous sections. The *Basic ACS* algorithm performs relatively well when there is some gradient information between different areas. In this image, the cyan and the blue contours (corresponding to the PZ and CZ, respectively) are practically attached to the borders; the contour corresponding to the TZ (yellow contour) is almost perfect. Nonetheless, the tumour area invades the TZ area, since the intensities contrast between them is not enough in order to define a good edge.

For the case of *Simplified ACS*, the results are not so satisfactory. This algorithm is uses simplifications of the formulas described in [12]. It has been tested for heart's left ventricle segmentation and it has a good performance in both, quality and speed, but for prostate segmentation, this is not true. Except for the TZ, all the results take part of the other zones, and only when the edge is too strong, the snake seems to be perfectly attached to it.

**Figure 12.** Segmentation results and ground truth for all the algorithms implemented.

Finally, as it is expected, the Gradient Vector Flow implementation overcomes all the problems, and the curves are strongly attached to the borders of each zone. There is only a little overlapping between the tumour and the TZ.

Another point to observe is the smoothness of the results from *Basic ACS* and *GVF ACS*. None of them include either a shape prior or the smoothing factor. This allows the snake to evolve in any shape, which is desired in this case, since there is not a single prior that work for all the prostate zones. Additionally, for the same reason, a generic post-processing that increases the smoothness of the curves is not trivial to implement. This problem is visible mainly in the result given by the *Basic ACS*, where all the contours do not present smooth borders. Nonetheless, the *GVF ACS* generates a result that even without these factors, the outcomes are not so jagged. Once again, the addition of the Gradient Vector Flow component surpasses the results given by the *Basic ACS*.

It is worth mentioning that the initialization (in all the cases) requires to be inside of the object to produce valid results. Since a lot of objects are present in an MR image, a larger contour than the object will produce poor segmentation results. In order to reduce this effect, the image to be segmented is cropped to an square of 150x150 pixels around the center of the initialization contour before starting the segmentation process. As a consequence, not only the influence of other objects over the snake is reduced, but also it reduces the execution time, since less operations are required at each iteration of the algorithm.

In order to measure quantitatively the efficiency of the three model implementations, a contour evaluation of this test was performed. For this evaluation, the Hausdorff distance (HD) and Dice index (DI) metrics for the various segmentation results were compared with a ground truth provided by a specialist. The Hausdorff distance between two finite set of points A and B is defined as follows:

$$HD(A, B) = \max(h(A, B), h(B, A))\tag{22}$$

where *h*(*A*, *B*) is the directed Hausdorff distance defined by:

$$h(A, B) = \max\_{a \in A} \min\_{b \in B} ||a - b||\tag{23}$$

where . is the Euclidean distance and *A* and *B* are the resulting contours of the segmentation.

In computer vision and image processing applications, the largest segmentation errors quantified using the Hausdorff Distance can be a good measure of the usefulness of the segmentation for the intended task. Intuitively, *HD* is the longest distance one has to travel from a point in one of the two sets to its closest point in the other set. In image segmentation, this metric is computed between boundaries of the estimated and the ground-truth segmentation, which consist of curves in 2D. We computed the Hausdorff distance between the generated segmentation and the ground truth. The methods with the lowest distance are considered to have yielded the best segmentation results.

Another popular metric for image segmentation tasks is based on the Dice coefficient, which is essentially a measure of overlapping between two samples and it is equivalent to the F1 score. This measure ranges from 0 to 1, where a Dice coefficient of 1 denotes perfect and complete overlap. This metric is calculated as:

$$\operatorname{Dice}(A, B) = \frac{2 \left| A \cap B \right|}{\left| A \cup B \right|} \tag{24}$$

In Table 1, we summarize the results obtained from the performance analysis for the three segmentation methods, using the two above-mentioned metrics. This procedure has been done for each of the prostate's gland zones, which corresponds to the examples of the annotations depicted on Figure 12. The ground truth is used as a reference image to calculate all the measurements. First, according to the initialization image, a comparison of the results between *Simplified ACS*, *Basic ACS* and *GVF ACS* is performed. It can easily be inferred from the table that the *Basic ACS* yields more accurate segmentation results in all the prostate zones than the *Simplified ACS*. Moreover, it can be seen that the *GVF ACS* implementation improves the results in all the cases.

**Table 1.** Analyses of Hausdorff distance (in mm) and Dice index for the CZ, PZ, TZ and tumour area (TUM).

The geometry of the *CZ* zone is considered the most difficult part of the prostate MR image because it is difficult sometimes to demarcate the edges of this area, given that it is surrounded by the other zones; for this reason, it is important to do a good initialization. The results for the *PZ* are good in all the cases, especially when GVF is used. Moreover, the results on the *TZ* area are good for all the algorithms considered in this study. Finally, the results for the tumor area indicate that the best segmentation is obtained using *GVF ACS* methods.

Finally, in order to assess the performance of the various studied segmentation methods, a time analysis has been performed. As in Table 1, the analysis has been divided by zones and by algorithm. Then, the execution time of the segmentation algorithm has been assessed for each zone, using the same initialization as in Figure 11. For this test, a computer with a microprocessor i7-8650U @ 1.90 GHz, 16 gb of RAM and running Linux 4.15.0-64 has been used. The execution times are depicted in Table 2. All the measurements has been done in milliseconds. Both *Simplified ACS* and *Basic ACS* run at maximum 200 iterations, and the *GVF ACS* runs at maximum 80 iterations. Both Tables 1 and 2 confirm what it has been said about the *Simplified ACS* algorithm. It is the least accurate, but it runs faster than the other ones. On the other hand, even though the *GVF ACS* runs less iterations than the *Basic ACS*, their running time are similar. This is due to the *GVF ACS* algorithm computes the same code, plus the additional ones in order to get the Gradient Vector Flow contribution. Nonetheless, for this work, this execution time could have been reduced if the equations related to the shape prior and the smoothness would have been removed. Only the weights have been set to zero for those priors. This code has not been erased so that other developers can prove other parameters values, and try

to find other combination of weights to make this program work with either, other image modalities, or other organs to be studied.

**Table 2.** Execution time (in ms) for each segmentation algorithm, and for each prostate zone.

These results show that a good segmentation can be obtained in a reasonably time, which is important for both clinical and machine learning applications. As it is expected, one result is not representative of the quality of these algorithms. An analysis of a larger population must be done in order to have a real measure of quantitative error of our implementations. The work we introduced here is to present the tool, the algorithms implemented and how to use them. As future work, we plan to perform a more comprehensive study encompassing other methods known to have good results to contrast results as well as other regions of interest within the prostate. This will enable us to assess the validity of the proposed segmentation approach in other setting, and to get feedback from a larger pool of specialists. Other avenues of research pertain the validation of the proposed tool, as well as of the segmentation algorithm in a larger pool o prostate MRI images, as well as with a wider number of clinicians. We also plan to integrate other segmentation methods into the proposed GUI tool to carry out a more comprehensive quantitative comparison and to investigate its diagnostic differentiation with other prostate diseases such as BPH, chronic or acute prostatitis.

#### **5. Conclusions**

In this paper we described the implementation of a snake-based segmentation algorithm, based on an Active Contours Models (ACM), which has been integrated into a GUI based application. We focused our research in a method originally conceived for segmenting the shape of the heart's left ventricle. A modification of this specific prior has yielded very good results in the context of prostate cancer segmentation: the snake models tries to evolve and to adjust to the boundaries of the proposed zones. The problem with this implementation is that the performance of this previous model highly depends on the quality of the input image, as it uses the pixel intensities and the edges of the object in order to determine the correct curve evolution. As it can be expected, if the boundaries are not well defined, or if the intensities of different sections have similar values, then either the model will not be able to assess the difference, preventing the algorithm of evolving, or it will cause the snake to overlap with other elements in the image. This scenario can be observed in Figure 12, where the red contour overlaps with the area that corresponds to the TZ.

Furthermore, in order to get that segmentation, several shapes for the initial contours had been tested to obtain the best results, which might be very different from one run of the algorithm to another. This entails that the implemented snake-based algorithms are very sensitive to the initialization conditions. This type of problems have fostered researches in image analysis and computer vision to find novel ways of constraining the shape of the force that pushes the objects' boundaries in the right direction. Doing so, the segmentation algorithm becomes robust against to initialization conditions and less sensitive to the quality of the image.

Regarding the performance of both implementations, the *Simplified ACS* and *Basic ACS*, their behaviour are not the same. Since *Basic ACS* includes all the weights and all the functions described by the original document, the boundaries of the final segmentation tends to be closer to the real boundaries of the object to segment than in the case of *Simplified ACS*. Nonetheless, both of them fail in the task of finding the boundaries when the image intensities between two contiguous areas are not noticeably different, like in the case of the tumor and the central zone in Figure 12. This problem has been partially overcame when using the Gradient Vector Flow external force.

Further improvements to this algorithm will thus be directed towards implementing other types of Gradient Vector Flow energy optimization algorithms, or exploring how to include shape priors for each zone of the prostate separately. These modifications would help to create a segmentation that does not go too far from where the real object is. Additionally, no pre-processing nor post-processing functions have been included. Adding pre-processing operations might help to remove the noise that addresses the snake towards the wrong directions, and post-processing of the resulting contour will help to make it smoother, and closer to the real boundaries of the object.

Moreover, since the DICOM files can contain different image modalities, and each modality has its particularities, a profiles section can be added to the program. Each profile will contain a set of parameters (priors weights, GVF weight, regularization constants, and iterations) specific for each algorithm. These profiles will be chosen in a way that they will produce the best segmentation results for the selected image modality, increasing the scope of application of the program.

Other avenues of research pertain the validation of the proposed tool, as well as of the segmentation algorithm in a larger pool o prostate MRI images, as well as with a wider number of clinicians. We also plan to integrate other segmentation methods into the proposed GUI tool to carry out a more comprehensive quantitative comparison and to investigate its diagnostic differentiation with other prostate diseases such as BPH, chronic or acute prostatitis

**Author Contributions:** Conceptualization, C.M. and J.R.; methodology, all authors.; software, J.R.; validation, C.M., G.O.-R.; resources , C.M and G.O.-R.; data curation, C.M.; writing—original draft preparation, J.R. and C.M.; writing—review and editing, all authors; visualization, all authors; supervision, C.M., G.O.-R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was partially funded by the Spanish R+D+I grant n. TIN2012-37171-C02-01, by UdG grant MPCUdG2016/022. The Regional Council of Burgundy under the PARI 1 scheme also sponsored this work. C. Mata held a Mediterranean Office for Youth mobility grant.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Medical Assistant Mobile Application for Diabetes Control by Simulating a Compartmental Model**

**Martín Hernández-Ordoñez 1,\*, Marco Aurelio Nuño-Maganda 2, Carlos Adrián Calles-Arriaga 2, Abelardo Rodríguez-León 3,\*, Guillermo Efren Ovando-Chacon 1,\*, Rolando Salazar-Hernández 4, Omar Montaño-Rivas 5 and José Margarito Canseco-Cortinas 2**

Received: 29 July 2020; Accepted: 25 September 2020; Published: 29 September 2020

#### **Featured Application: Mobile application for iterative simulation of a diabetic model.**

**Abstract:** This paper presents an educational mobile assistant application for type 1 diabetes patients. The proposed application is based on four mathematical models that describe the glucose-insulin-glucagon dynamics using a compartmental model, with additional equations to reproduce aerobic exercise, gastric glucose absorption by the gut, and subcutaneous insulin absorption. The medical assistant was implemented in Java and deployed and validated on several smartphones with Android OS. Multiple daily doses can be simulated to perform intensive insulin therapy. As a result, the proposed application shows the influence of exercise periods, food intakes, and insulin treatments on the glucose concentrations. Four parameter variations are studied, and their corresponding glucose concentration plots are obtained, which show agreement with simulators of the state of the art. The developed application is focused on type-1 diabetes, but this can be extended to consider type-2 diabetes by modifying the current mathematical models.

**Keywords:** mathematical models; Iterative simulation; compartmental model; diabetes control; mobile assistant

#### **1. Introduction**

In recent years, obesity and overweight have become a severe problem for human health. In 2014 more than 1.9 billion adults aged 18 and older were overweight globally, from which over 600 million were obese. The data indicated that 39% of these adults were overweight and 13% were obese [1]. Body mass index (BMI) provides the most useful population-level measure of overweight and obesity as it is the same for both sexes and all ages of adults. Patients with higher BMI are at greater risk for having diabetes. In 2014, more than 422 million people had diabetes in the world, with a global prevalence of 8.5% among adults over 18 years old. Diabetes is a leading cause of blindness, kidney failure, heart attacks, stroke and lower limb amputation. In 2012, an estimated 1.5 million deaths were directly caused by diabetes, and another 2.2 million deaths were attributable to high blood

glucose [2]. Diabetes was promoted from 21st place in 1990 to the 14th in 2010 in the list of death causes worldwide [3], and recently, to the 6th place in 2016 [4].

In Mexico, about 72.5% of adults from 18 years and over, 36.3% of teenagers, and 33.2% of children were overweight or obese in 2016. Diabetes statistics are not too optimistic; Mexico is in sixth place in countries with more diabetes, just behind China, India, US, Brazil and Russia [5]. Recent statistics show that 9.17% of adults from 18 years (about 6.4 million people) have diabetes [6]. Type 2 diabetes has become the leading cause of death in Mexico. The disease claims nearly 80,000 lives each year, and forecasters say the health problem is expected to get worse in the decades to come [7].

On the other hand, smartphones have been accessible to more people due to the high worldwide market growth rate. In 2016, around 1.5 billion smartphones were sold worldwide [8]; by the end of 2016, more than 3.2 billion individuals worldwide used the Internet regularly, accounting for nearly 45% of the world's population; about 2.5 billion people access the web via mobile phone [9], which represents novel opportunities for designing applications running on smartphones, accessible to large population groups. Specifically, healthcare and disease management have benefited with applications such as fitness tracking, drug dose calculators, information apps, tutorials, surveys, diaries, and other tools created to help empower users to put more emphasis on their health care. Some benefits of using smartphones with Internet connectivity consists of the possibility of sharing data with family, friends, classmates, colleges, collaborators but especially with professionals who can guide the user and give feedback regarding treatment to be followed.

This paper is structured as follows. In Section 2, a concise review of the state of the art is given. A description of the models used in the developed application is provided in Section 3. In Section 4, the results and discussion are described. Finally, in Section 5, the conclusion is presented where final recalls are established, and future work is outlined.

#### **2. State of the Art**