Francisco Durán

Camilo Rocha

Gwen Salaün

A Rewriting Logic Approach to Resource Allocation Analysis in Business Process Models

Keywords: Business processes, BPMN, resource allocation, rewriting logic, executable specification, automated verification, Maude

This paper presents an approach for the modeling and analysis of resource allocation for business processes. It enables the automatic computation of measures for precisely identifying and optimizing the allocation of resources in business processes, including resource usage over time. The proposed analysis, especially suited to support decision-making strategies, is illustrated with a case study of a parcel ordering and delivery by drones that is developed throughout the paper. The paper comprises an encoding of a significant and expressive subset of the Business Process Model and Notation (BPMN) in rewriting logic, an executable logic of concurrent change that can naturally deal with state and with concurrent computations. The encoding is by itself a formal semantics and interpreter of the BPMN subset that captures all concurrent behavior and thus is used to simulate the concurrent evolution of any business process with a given number of resources and replicas.

Introduction

Business process optimization is a strategic activity in organizations because of its potential to increase profit margins and reduce operational costs. One of the main challenges in this activity is concerned with the problem of optimizing allocation and sharing of resources. This is a key question because having a clear answer on how resources are used can help in identifying availability, detect bottlenecks, and improve sharing in order to make processes more efficient. However, analyzing resource usage in business processes is far from being an easy task. First, there is the need for a formal model expressive enough to capture the semantics of processes. Second, resources have to be included in the formal model supporting multiple concurrent executions of a process. Last, automatic verification is required to compute certain metrics and properties of interest (e.g., resource charge, resource occupancy, or usage percentage of each resource replica) on hundreds and even thousands of simulations of a process with its resources.

This paper presents a formal executable model for business processes with resource allocation and automatic verification techniques for analyzing resource usage. The formal model relies on an encoding of a subset of the Business Process Model and Notation (BPMN) in rewriting logic [START_REF] Meseguer | Conditional Rewriting Logic as a Unified Model of Concurrency[END_REF]. This specification is executable in Maude [START_REF] Clavel | All About Maude -A High-Performance Logical Framework, How to Specify, Program and Verify Systems in Rewriting Logic[END_REF], and supports the concurrent execution of a process with different types of resources and with multiple replicas on any given workload. The verification techniques use Maude's rewriting tools for evaluating expected values of any numerical or path expression in the executable model. The Maude specification supports BPMN activity and collaboration diagrams, four types of gateways, quantitative aspects (such as times associated to flows and tasks, as well as probabilities associated to branching behaviour), loops, resource description, and unbalanced processes (which are workflows without a strict correspondence between split and merge gateways).

The overall idea, from a verification viewpoint, is that multiple concurrent executions of a process compete for the shared resources. With the help of Maude, this approach enables the automatic analysis of how resource usage evolves over time when varying the workload and the amount of resources. This analysis is performed on design models, without the need of an implementation of the system running on real resources. All that is needed is the BPMN specification of the process workflow completed with expert information on the operation of the system. The usefulness of this kind of analysis is illustrated with several experiments on real-world processes, confirming that the approach can help in detecting resource usage problems, thus ultimately leading to the improvement of the business process by optimizing its resource allocation. In particular, the experiments presented in this paper identify, e.g., low-level occupancy of resources and undesirable patterns of resource usage, such as sequential dependencies and bottlenecks provoked by some highly used resources that may induce performance fall-downs. Furthermore, a recommender system is proposed to find an optimal assignment of resources following Monte Carlo simulations of the business process.

In summary, the main contributions of this work are: (i) a formal semantics captured by an encoding in rewriting logic of a significant subset of BPMN, including activity and collaboration diagrams, four types of gateways, unbalanced workflows, loops, resource description, and concurrent process execution; (ii) automated formal analysis of resource allocation properties (e.g., accumulated time of occupancy for resources) performed at design time and a recommender system computing the optimal assignment of resources; (iii) the evaluation of the approach on a workbench of real-world BPMN processes. An early version of this paper has been published in [START_REF] Durán | Computing the parallelism degree of timed BPMN processes[END_REF] and has been significantly extended here as follows:

(i) a version of BPMN with explicit description of resources is considered;

(ii) the Maude encoding has been extended to support resources and the multiple execution of a process (tokens, workload, scheduling) required for simulating the evolution of resource usage over time; (iii) several properties of interest regarding average execution time and resource occupancy, that can be analysed using the proposed framework, have been identified; (iv) the approach has been applied to several case studies for validation purposes; and (v) the introduction, related work section, and conclusion have been fully revised, by taking the aforementioned extensions and improvements into account. The rest of the paper is organized as follows. Section 2 introduces the main features of BPMN. Section 3 briefly presents rewriting logic and explains the encoding of the considered subset of BPMN. Section 4 presents how resource analysis is achieved using Maude's rewriting-based tools. This section also describes several experiments and illustrates how the automatic verification can help in improving process definitions. Section 5 surveys related work and Section 6 concludes the paper.

Business Process Model and Notation

BPMN 2.0 (BPMN, as a shorthand, in the rest of this manuscript) was published as an ISO/IEC standard [START_REF]International Standard 19510, Information technology -Business Process Model and Notation[END_REF] in 2013 and is nowadays extensively used for modeling business processes. In this paper, two main kinds of BPMN diagrams are considered, namely, activity and collaboration diagrams. In particular, the focus is on the BPMN elements related to control-flow modeling and behavioral aspects that can be represented in BPMN constructs. Beyond those constructs, resource description and allocation are also considered, for which an automated analysis method is proposed.

Figure 1 summarizes the BPMN constructs supported in this work. These elements are used to develop activity and collaboration diagrams of process models. In addition to the description of specific tasks and their sequencing, collaboration diagrams also involve pools and lanes, which are structuring elements that split processes into pieces.

Specifically, the node types event, task, and gateway, and the edge type sequence flow are considered. Start and end events are used, respectively, to initialize and terminate processes. A task represents an atomic activity that has exactly one incoming and one outgoing flow. A task may have a duration (expressed as a stochastic expression) and may produce an event message. A sequence flow describes two nodes executed one after the other, i.e., by imposing an execution order with possible delays. The timing information associated to tasks and flows (durations or delays) is described either as a literal value (a non-negative real number, possibly 0) or sampled from a probability distribution function according to some meaningful parameters. The probability distribution functions currently available include exponential, normal/Gauss, and uniform (see, e.g., [START_REF] Walck | Hand-Book on Statistical Distributions for Experimentalists[END_REF]).

Gateways are used to control the divergence and convergence of the execution flow. Four types of gateways are considered: exclusive, inclusive, parallel, and event-based. Gateways with one incoming branch and multiple outgoing branches are called splits, e.g., split inclusive gateway. Gateways with one outgoing branch and multiple incoming branches are called merges, e.g., merge parallel gateway. An exclusive gateway chooses one out of a set of mutually exclusive alternative incoming or outgoing branches. For an inclusive gateway, any positive number of branches among all its incoming or outgoing branches may be taken (both BPMN 1.0 and 2.0 semantics for inclusive gateways are supported). A parallel gateway creates concurrent flows for all its outgoing branches or synchronizes concurrent flows for all its incoming branches; event-based gateway takes one of its outgoing branches or accepts one of its incoming branches based on events. Event-based split gateways may have a default branch fired by a timeout.

Data-based conditions for split gateways are modeled using probabilities associated to outgoing flows of exclusive and inclusive split gateways. The probabilities of the outgoing flows in an exclusive split must sum up to 1, while each outgoing flow in an inclusive split can be equipped with a probability between 0 and 1 without a restriction on their total sum. Workflows with looping behavior are supported, as well as unbalanced workflows.

Each lane in a collaboration diagram corresponds to a specific role or to a specific resource. Collaboration diagrams are also often used for modeling distributed systems and, in that sense, they heavily rely on message events and event-based gateways. Once completed, a task may generate a message with a flow in another lane as target. An event-based split gateway is triggered whenever a message is available from one of its outgoing flows or when its timer event has completed.

Instead of implicitly associating resources to lanes, resources are explicitly defined at the task level. A task that requires resources can include, as part of its specification, the number of required instances (or replicas) of a resource. Indeed, several tasks could compete for the same resources. Furthermore, since multiple instances of a same process may be executed concurrently, different instances may also access and thus compete for the shared resources. Section 3 presents how the number of available resources can be specified and Section 4 how the optimal number of resources can be computed.

The execution semantics of BPMN is informally described in official documents such as [START_REF] Omg | Business process model and notation[END_REF][START_REF]International Standard 19510, Information technology -Business Process Model and Notation[END_REF]. Although the focus here is on the provided analysis techniques, this work also gives a formal semantics to BPMN with time and resources. Last but not least, in this paper, BPMN processes are assumed to be syntactically correct. This can be enforced using existing works and tools, e.g., [START_REF] El-Saber | BPMN formalization and verification using maude[END_REF], the Activiti BPM platform, Bonita BPM, or the Eclipse BPMN Designer.

Running example. To introduce and illustrate the use of these elements and the approach presented in this paper, Figure 2 presents a process describing a parcel ordering and delivery by drones. This BPMN process is presented as a collaboration diagram consisting of three lanes, one for the client, one for the order management and one for the delivery process. In this example, the client first signs in and then repeatedly looks for products. Eventually, the client can decide to give up (termination) or to make an order by submitting it to the order management lane. The client then waits for a response (acceptance or refusal of this order). If the order can be completed, the client pays for it and then receives the parcel. Otherwise (timeout or order refused), the client fills in a feedback form. As far as the management lane is concerned, the first task aims at verifying whether the goods ordered by the client are available. If they are not available, then the order is canceled; otherwise, the order is confirmed. The order management takes care of the payment of the order whereas the delivery lane is triggered to prepare the parcel to be delivered by a drone. This process exhibits different kinds of gateways, probabilities for choice gateways, stochastic functions for time associated to tasks, a loop (Search products task), and unbalanced structures (for the eventbased gateways used in the client lane and the exclusive gateway after the initial availability check).

To simplify the exposition of the running example (Figure 2), the delays in all flows are set to 0 and the specification of the task duration has been placed at the bottom-left corner of the process description. For instance, the duration of the Sign in task follows a normal distribution with mean 1 and variance 0.5, and the Search products task follows a uniform distribution in the interval [START_REF] Arévalo | A metamodel to integrate business processes time perspective in BPMN 2.0[END_REF][START_REF] Poizat | Checking business process evolution[END_REF]. Exclusive split gateways are modeled using probabilities associated to outgoing flows. For instance, notice the exclusive split after the Search products task in the Client lane of the running example, which has outgoing branches with probabilities 0.6, 0.2, and 0.2, specifying the likelihood of following each corresponding path. All these numbers (durations and probabilities) can be provided by experts or can be obtained by executing the process in practice on a reasonable period of time (a few days) and by extracting this information from the generated logs.

The process in Figure 2 relies on employees for parcel packing and drones for parcel delivery. Notice the small circles at the bottom-right corner of the Prepare parcel and Deliver parcel tasks, indicating that one instance of the employee resource and another one of the drone resource are required, respectively, for the tasks completion.

From BPMN to Maude

This section presents the encoding of BPMN processes including time and resources in Maude. The Maude system is chosen as the formal specification language and verification platform for several reasons. Its language is expressive enough for specifying all the aspects of BPMN introduced in the former section. It also offers several rewriting-based tools that are adequate for simulating the specification, and for computing numerical results for timing properties and resource analysis. This section first gives an overview of rewriting logic and the Maude framework. The Maude specification of BPMN is then described in two parts: the encoding of the process structure and the description of the semantics of the BPMN subset using rewrite rules. The complete Maude implementation is available online at [START_REF] Durán | A Note on Resource Allocation Analysis of BPMN Processes[END_REF].

Rewriting Logic and Maude Overview

Rewriting logic [START_REF] Meseguer | Conditional Rewriting Logic as a Unified Model of Concurrency[END_REF] is a logic of change that can naturally deal with state and with highly nondeterministic concurrent computations. A rewrite logic theory is a tuple (Σ, E ∪ A, R), where (Σ, E ∪ A) is a membership equational logic [START_REF] Bouhoula | Specification and proof in membership equational logic[END_REF] theory with Σ its signature, E a set of conditional equations and sort membership axioms, A a set of equational axioms (e.g., associativity, commutativity and identity) so that rewriting is performed modulo A, and R is a set of labeled conditional rules.

Maude [START_REF] Clavel | All About Maude -A High-Performance Logical Framework, How to Specify, Program and Verify Systems in Rewriting Logic[END_REF] is a high-level language and a high-performance interpreter that supports membership equational logic and rewriting logic specification and programming of systems. Thus, Maude integrates an equational style of functional programming with rewriting logic computation. Thanks to its efficient rewriting engine and its metalanguage capabilities, Maude turns out to be an excellent tool for creating executable environments of various logics, models of computation, theorem provers, or even programming languages.

A functional specification must be terminating, confluent, and sort-decreasing. Computation in a functional module is accomplished by using the equations as simplification rules from left to right until a canonical form is found. Some equations, such as those expressing the commutativity of binary operators, are -however-not terminating. Nonetheless, they are supported by means of operator attributes, so that Maude performs simplification modulo the equational theories provided by such attributes, which can be associativity (assoc), commutativity (comm), identity (id), and idempotency (idem). The above properties must therefore be understood in the more general context of simplification modulo such equational theories.

In Maude, a distributed system is axiomatized by a rewrite theory describing its states as an algebraic data type (an equational sub-specification) and a collection of conditional rewrite rules specifying its behavior. Rewrite rules are written crl [l] : t => t if C, with l the rule label, t and t terms, and C a guard or condition. Rewrite specifications are not required to be terminating nor confluent. Rules describe the local, concurrent transitions that are possible in the system, i.e., when a part of the system state fits the pattern t, then it can be replaced by the corresponding instantiation of t . The guard C acts as a blocking precondition: a conditional rule can only be fired if its condition is satisfied. Sometimes rules are given without label or condition (which can be assumed to be true).

In the Maude language, object-oriented systems can be specified by objectoriented modules in which classes and subclasses are declared, with the usual support for inheritance, dynamic binding, etc. A class is declared with syntax class C | a 1 : S 1 , . . . , a n : S n , where C is the name of the class, a i are attribute identifiers, and S i are the sorts of the corresponding attributes. The objects of a class C are then record-like structures of the form < O : C | a 1 : v 1 , . . . , a n : v n >, where O is the name of the object and v i are the current values of its attributes.

In a concurrent object-oriented system, the concurrent state, which is called a configuration, consist of a multiset of objects and messages. Rewrite rules then define transitions between such configurations. These transitions represent the different actions that may occur in the system. For instance, there will be rules modeling the effects of events, or the synchronous and asynchronous communication events of objects and messages. The general form of a rewrite rule r is the following:

crl [r] : < O 1 : C 1 | atts 1 > ... < O n : C n | atts n > M 1 ... M m => < O i 1 : C i 1 | atts i 1 > ... < O i k : C i k | atts i k > < Q 1 : C 1 | atts 1 > ... < Q p : C p | atts p > M 1 ... M q if Cond .

where

• r is the rule label,

• M 1 ...M m and M 1 ...M q are messages,

• O 1 ...O n and Q 1 ...Q p are object identifiers, • C 1 ...C n , C i 1 ...C i k and C 1 ...C p are classes, i 1 ...i k is a subset of 1...

n, and

• Cond is a Boolean condition (the rule's guard). The result of applying such a rule is that:

• messages M 1 ...M m disappear, i.e., they are consumed,

• the state, and possibly the classes of objects O i 1 ...O i k may change,

• all the other objects O j vanish,

• new objects Q 1 ...Q p are created, and

• new messages M 1 ...M q are created, i.e., they are sent. The real-time aspects are modeled using Real-Time Maude [START_REF] Ölveczky | Semantics and pragmatics of real-time maude[END_REF], which supports the formal specification and analysis of real-time systems. Specifically, Real-Time Maude provides a sort Time to model the time domain, which in this work is assumed to be dense (i.e., represented as rational numbers). Then, given a system configuration and a time amount, time elapse is modeled with tick rules like

crl [l] : { t, T } => { t , T + τ } if C .
where t and t are system states, T is the global time, and τ is a term of sort Time that denotes the duration of the rewrite, and that affects the global time elapse. Since tick rules affect the global time, in Real-Time Maude time elapse is usually modeled by one single tick rule and the system dynamic behavior by instantaneous transitions [START_REF] Ölveczky | Semantics and pragmatics of real-time maude[END_REF]. Although there are other sampling strategies, in the most convenient one this single tick rule models time elapse by using two functions: the delta function, that defines the effect of time elapse over every model element; and the mte (maximal time elapse) function, that defines the maximum amount of time that can elapse before any action is performed. Then, time advances non-deterministically by any amount τ, which must be less than or equal to the maximum time elapse of the system.

Note that the above is a general form for the tick rules. As it will be seen in the following sections, the concrete representation of states and global time depend on the concrete system specification. This simple approach allows for having multiple clocks and timers, and model the different real-time events that could happen in the system quite naturally.

Process Description

In the Maude specification of BPMN, a process is represented as an object with sets of flows and nodes as attributes. Nodes can be of five different types: start, end, task, split, or merge. The representation of each of these types of elements includes the necessary information. A task node involves an identifier, a description, two flow identifiers (input and output), a stochastic function modeling its duration (0 if there is no duration), a set of resources required for its execution, and a set of messages to be delivered after its completion. A split node includes a node identifier, a gateway type (exclusive, parallel, inclusive, or event-based), an input flow identifier, and a set of output flow identifiers. A merge node includes a node identifier, a gateway type, a set of input flow identifiers, and an output flow identifier. The representation of a flow includes a probability distribution function specifying its delay, a message produced by a task that blocks the flow until the message is received, and a timer representing a delay after which the execution can be triggered.

Figure 3 gives an excerpt of the representation for the running example. It shows how a Process object has attributes with the definition of its nodes and flows connecting them. For example, the exclusive split g2 has as incoming flow cf4 and outgoing flows cf5, cf6, and cf7, with associated probabilities 0.6, 0.2, and 0.2, respectively. As another example, the event-based split gate g3 has as < pid : Process | nodes : (start(initial, cf1), merge(g1, exclusive, (cf2, cf5), cf3), split(g2, exclusive, cf4, ((cf5, 0.6) (cf6, 0.2) (cf7, 0.2))), split(g3, eventbased, cf8, (cf9, cf10, cf11)), task(t10, "Prepare parcel", mf7, df1, Norm(5.0, 4.0), employee, empty), task(t11, "Deliver parcel", df1, df2, Unif(5.0, 30.0), drone, parceldelivered), ...), flows : (flow(cf1, 0), flow(cf9, 0, message(orderconfirmed, "Order confirmed")), flow(cf10, 0, message(ordercanceled, "Order canceled")), flow(cf11, 0, timer(timeout, 60)), ...) > incoming flow cf8 and outgoing flows cf9, cf10, and cf11. These flows are defined in the set of flows.

The transformation from the BPMN diagrammatic representation of processes into the corresponding Maude representation is automated by an extension of the VBPMN platform [START_REF] Krishna | VBPMN: automated verification of BPMN processes (tool paper)[END_REF].

Execution semantics

The operational semantics of BPMN is defined using rewrite rules, modeling how tokens evolve through a process, thus defining the execution semantics of BPMN. Each observable action is modeled as a rewrite rule. For instance, when a token arrives at a parallel split gateway, the token corresponding to the incoming flow is removed, and one token is added for each outgoing flow. Technically, rewrite rules operate on systems composed of a process object and a simulation object.

Simulation object. While the process object introduced in Section 3.2 represents the BPMN process and does not change during an execution, the simulation object keeps information on the execution of the process. It stores a collection of tokens (in a scheduler, see below), a global time (gtime), a set of events (messages and timers), and a set of resources. It also keeps track of the quantities being measured during the analysis of a process. Figure 4 presents the structure of the Simulation object.

Tokens. Tokens are used to represent the evolution of the workflow under execution. A token is represented as a term token(TId, Id, T). Since several executions are simultaneously happening, each execution has a unique identifier, and tokens < s : Simulation | tokens : ..., ----scheduler gtime : ..., ----global time resources : ..., ----resource set events : ..., ----event set process-execs : ...., ----execution times sync-times : ..., ----synchronization times task-times : ..., ----task execution times ... > are identified by the execution instance TId they belong to, and the flow or node Id they are attached to. The expression T represents a timer, of sort Time, modeling a delay on the token. Once this timer becomes 0, the token may be consumed.

Scheduling. Tokens are stored in a scheduler implemented as a priority queue, so that they are kept according to their due time. However, even with its timer set to 0, the token at the front may be not enough to fire some action. Consider, for example, a task that requires some resource that is not available or a parallel merge for which some incoming flow is not yet active. To avoid blocking situations, the scheduler is provided with a shifting mechanism, which moves the first active token to the front of the scheduler in case the current head cannot fire the corresponding action. This scheduler is similar to those used in typical discrete event simulations.

Events. A message event may be associated to a flow, and the flow is blocked until the message is received. A timer event may be associated to a flow and has a delay as parameter. When the flow of execution arrives at a timer event, its countdown is started: once the countdown is completed, the token moves to the outgoing flow. Both message and timer events are usually associated to eventbased gateways, but it is not necessarily the case (see, e.g., the initial flow for the order management lane in the process in Figure 2). Asynchronous events are modeled using an event set in the Simulation object. When a message is dispatched, a corresponding event is added to the set. Flows and gateways that are waiting for specific messages use this set to check whether the messages have arrived.

Resources. Each resource is represented with an identifier, the number of available replicas (initially the total number), the total amount of time this resource has been in use, and the intervals of time on which it was used. These two last parameters are required for analysis purposes only. When a task requires several resources, it atomically uses all of them at once, or waits for them to become available. Workloads. Simulation-based analysis techniques are typically parameterized by the workload that represents the way a system is used. They define the rate at which new instances of a given process are executed. Currently, closed workloads can be handled by specifying the number of executions and the rate at which executions are started, that is, their inter-arrival times. Both the number and the rate are specified as stochastic expressions.

Rewrite rules for BPMN constructs. Rewriting rules represent how tokens evolve through the process and events are fired, thus defining the execution semantics of BPMN. Each action supported by the system is modeled as a rewrite rule. These rules are overviewed in the rest of this section to gather an intuition on the formal semantics (see [START_REF] Durán | A Note on Resource Allocation Analysis of BPMN Processes[END_REF] for the complete specification).

Start/end events. Figure 5 depicts the rule for the start event. When there is a token in the execution TId in the start node NId with delay 0 (note the token at the front of the scheduler in the Simulation object in line 5), then this rule generates a new token on the outgoing flow of the selected node to initiate the execution of a process instance (line 10). The insert function puts this token in the scheduler and the eval function evaluates the stochastic expression SE specifying the delay of the outgoing flow FId to be assigned to the new token. Details on the initialization of time stamps and recorded times for the initiated execution have been replaced by ellipses. A termination rule, associated to stop events, consumes tokens when they arrive at those events.

Tasks. A task execution is modeled with two rules. The first rule, the init-Task rule shown in Figure 6, represents the task initiation, which is applied when a token with zero time is available for the incoming flow (line 5). If all the re- sources required by this task are available, which is checked with the allResource-sAvailable function (line 8), then a new token is generated with the task identifier and the task duration (line 12). Otherwise, the scheduler's token shifting mechanism is invoked (line 20). If available, all required resources are removed from the set of resources, and the time those resources have been in use is updated (grabResources&updateTime function, line 18). Since all auxiliary functions in the right-hand side of the initTask rule are defined equationally, the checking and grabbing of resources are performed atomically, without introducing any blocking issues. Note also that rules update the information on execution times, task durations, etc. (see, e.g., the update of the task-tstamps attribute, lines [START_REF] Durán | Stochastic analysis of BPMN with time in rewriting logic[END_REF][START_REF] Durán | Verifying timed BPMN processes using maude[END_REF][START_REF] El-Saber | BPMN formalization and verification using maude[END_REF][START_REF] Gagné | Time-bpmn[END_REF]). This information is important for analysis purposes, as it will be seen in Section 4.

A second rule, which models task completion, is triggered when there is a token for that task with zero time. In that case, the token is consumed and a new one is generated for the outgoing flow. All resources are released, and all the message events associated to that task, if any, are added to the events set.

Exclusive gateways. There are two rules for the exclusive gateways, namely, one for the split and one for the merge. The rule for the split applies when a token with zero time is available on its incoming flow. A uniformly sampled probability distribution is used to choose the branch to be executed. The newly created token is assigned with its run-to-completion time generated by evaluating the stochastic expression associated to the chosen outgoing flow-this is actually the case every time a new token is added for a flow. The exclusive merge gateway is triggered when one of its incoming flows has a token with zero time. In that case, a new token is generated, assigned to the outgoing flow, and added to the scheduler.

Parallel gateways. The parallel split gateway rule is triggered when a token with zero time corresponding to the input flow is available. If so, the token is consumed and one token is added to each of its outgoing flows. The merge rule for the parallel gateway is executed when there is a token with zero time for each incoming branch. In that case, these tokens are removed and a new token is generated for the outgoing flow. In the merge rule, synchronization times are also updated.

Inclusive gateways. The split rule applies when a token with zero time is available at the incoming flow. Since all outgoing branches are equipped with probabilities, a function in charge of computing the subset of branches to be triggered is invoked. For each one of the selected branches, a new token is added to the scheduler. Regarding merge gateways, both BPMN 1.0 and 2.0 semantics are supported in this research. In BPMN 2.0, merge inclusive gateways behave like exclusive ones. The 1.0 version of the semantics is more involved [START_REF] Christiansen | Formal semantics and implementation of BPMN 2.0 inclusive gateways[END_REF], since the merge rule for the inclusive gateway is executed when all the expected tokens are available with zero time. This requires a global analysis. To check whether all expected tokens have arrived, a backwards traversal that explores the process upstream and checks whether there are tokens on their way to that merge is performed. In both cases, once the merge gateway is triggered, the incoming tokens are removed, a new token is added to the scheduler for the outgoing flow, and simulation information is updated with synchronization times.

Event-based gateways. When a token arrives at an event-based split gateway, the token is made active with its optional timer. In that rule, if there is an outgoing flow with a timer, an event is added with the corresponding time to the set of available events. Two additional rules handle the possible cases. If there is an outgoing flow with a message in the set of events, then that branch can be activated and one token is added for that flow. This behavior is formalized with the splitGatewayEventBased rule shown in Figure 7, where an event msg(MId) is in the events set and MId is the message associated to the flow FId1, which is one of the outgoing flows of the split NId. A second rule handles the case in which none of the expected messages is in the set of events, but there is at least one timer for which the time has expired. In this situation, a token is added to the outgoing flow corresponding to that timer. Loops and unbalanced workflows. The modeling of the BPMN execution semantics using tokens and their circulation through the process structure supports intricate constructs such as loops and unbalanced workflows. As far as looping behavior is concerned, a token may circulate back to an already visited flow without any additional treatment. Similarly, tokens can advance through flows that are part of balanced or unbalanced gateways, independently of their structure.

Resource Allocation Analysis

This section illustrates how resource allocation analysis can be performed with the proposed approach using the running example. It also includes a discussion on how the analysis can help in refining and improving resource allocation. Specifically, given the process description, a specification of resources and a workload, the experiments illustrate how information on execution times and resource usage is collected. This information is used to find the optimal allocation of resources that minimizes costs and execution times relative to an optimization goal. The interested reader is referred to [START_REF] Durán | A Note on Resource Allocation Analysis of BPMN Processes[END_REF] for further details on the experiments carried out on the running example.

Properties

The BPMN subset encoded in Maude is quite expressive and several kinds of properties can be computed, including timing and resource-based properties. These properties are meaningful when executing multiple instances of a process that compete for the shared resources. As for timing properties, the approach presented in this paper allows the computation of average execution times (AET) of a process, its variance (Var), and the average synchronization time (AST) for merge gateways, representing the time elapse from the arrival of the first token through one of its incoming flows to its activation. Synchronization times make sense only for parallel and BPMN 1.0 inclusive gateways, since there is no waiting nor synchronization time for the other gateways.

As far as resource-based properties are concerned, which is the main focus in this work, the following properties are computed:

• The global time usage of all instances of each resource R (GTU R ). E.g., when executing 10 instances of a process P, with an AET of 42, it is possible that the two instances of a resource A are used for 56 time units and the three instances of resource B for 60 time units. • The expression GTU 1 R denotes the average GTU of resource R (i.e., the GTU per instance of resource R). Thus, although in the previous example GTU B is greater than GTU A , GTU 1 A is 28 and GTU 1 B is 20. • The average usage percentage UP R for a resource R over the global execution time. E.g., continuing with the running example, in average, an instance of the resource A is used 24% of the global execution time when executing 200 instances of a process P. To verify these properties, Maude rewriting capabilities are used in order to simulate and extract analysis results on a given BPMN process. The simulation object presented in Section 3 is used to accumulate information of synchronization times, task duration, and resource usage. At the end of all executions, these results are used for computing the expected average times and resource usage percentages. Since the analyzed processes are assumed syntactically correct and processes that may lead to non-terminating analysis are not considered (e.g., loops without end events), the verification process always terminates. This is the case because all splits are probabilistic, and time duration and probabilities assigned to the branches respect specific assumptions (e.g., all probabilities are between 0 and 1, they sum up to 1 in exclusive branches, and times are positive).

Evaluation

Table 1 summarizes experimental results on execution times and resource usage on the parcel order and delivery example (Section 2). They were carried out on an iMac with 3,2GHz Intel Core i5 and 8GB of RAM. All simulations were performed assuming a given workload with a number of instances (first column) and an exponentially distributed interarrival time (λ = 4). Columns 2 to 6 contain, respectively, the average execution time (AET), its variance (Var), the average synchronization time for the parallel merge at the end of the delivery process lane (AST g8 ), the average synchronization time for the end events (AST ee ), and the total time to complete the execution of all instances. The next six columns show results on resource usage for employees and drones-note the use of subindices e and d for employees and drones, respectively. The final column gives the overall time needed to complete the analysis. All times are logical (time units), except for the ones in the last column that are given in seconds. Other information, such as the duration of each task and the synchronization time of each merge gateway is also collected.

These experiments consist of 100, 200, 400, 800, and 1600 instances for 2 employees and 3 drones. Note that the average execution and synchronization times clearly increase with the number of instances. This is because the more tokens compete for resources, the more time it takes to execute the process and for the tokens to reach the synchronization points. Note the relationship between AET and AST ee times, showing an unbalance between the three lanes: the client lane terminates earlier than the two other lanes, which exhibit a bottleneck because of the demand on the resources.

The global time (GTU) for each or all the instances increases with the number of executed instances. These times are particularly interesting because they can be materialized as costs (e.g., cost of a resource, salary of an employee/all employees). In relation with usage percentage (UP), the results indicate that the employees are "underused" since they work around 40% of the time, in contrast to the drones that are constantly busy and used about 90% of the time for delivering parcels. This may suggest an inappropriate allocation of resources. It is worth observing that, although the number of instances clearly affects all computed times, the results for resource usage (UP) are quite stable and a small number of instances is enough for obtaining a good approximation of these percentages.

Resource allocation impacts execution times (AET) and resource usage (UP) of a process. Figure 8 focuses on average execution time and depicts the results when the number of employees and drones vary for a fixed number of executions (400). The objective here is to reduce the average execution time for completing the process: the quicker the parcel is delivered, the more satisfied the client is. It can be observed that, independently of the number of employees, execution times are not satisfactory with 1 or 2 drones (between 400 and 800 time units). The time becomes reasonable for more than 3 drones (less than 300 time units) and tends to stabilize. It is also worth noting that, given its low usage rate, the number of employees does not impact significantly the execution time. For more than six drones, only going from one to two employees makes a significant impact in the AET values. Figure 9 gives a different point of view of resource usage by concentrating on each resource instance. Figure 9 (left) shows that employees are close to 100% usage only if there is 1 or 2 instances of that resource and at least 4 or 8 drones, respectively. If the number of employees increases, the usage percentage quickly drops, reaching a low level (e.g., 14% for 4 employees and 2 drones). This per-centage slightly increases with the number of drone instances (e.g., 34% for 4 employees and 5 drones), but remains low (around 30%). Figure 9 (right) shows that the drone usage is always quite high whatever the number of employees is. With only 1 or 2 drones, the usage percentage is almost at 100% and slightly decreases with 4 drones. When there are 6 drones and 1 employee, the percentage is still about 60%. Another interesting fact is that the number of employees barely impacts the drone usage percentage. For example, with 4 drones, the usage percentage is around 90% for any number of employees (varying from 1 to 10).

Optimal Resource Allocation

The collected data on execution times and resource usage can be used to compute the optimal allocation of resources for the running example. Given the data represented in Figures 8 and9, an optimization problem can be expressed as follows. Given expressions f cost and f aet representing the normalized functions for drone and employee cost (per time unit) and for average execution time, and given w i with i ∈ {cost, aet} their respective weights, the following multi-objective optimization problem min

x∈X i∈{cost, aet} w i f i (x),
for which Pareto optimal solutions are required, can be solved. For instance, assuming that the cost per hour of a drone and an employee is 20e and 50e, respectively, Figure 10 Table 2 shows the results of the simulations for the running example, with 3 employees and 8 drones. If the results in Tables 1 and2 are compared, an improvement in all aspects can be observed, as expected. Specifically, average execution times keep variance values and synchronization times at lower rates. With this new assignment of resources, AETs have dropped significantly. E.g., from 910 to 302 for the 1600-instances case. The total time has similarly been reduced in all cases, e.g., from 4263 down to 1627 in the 1600-instances case. To conclude, this improvement is a consequence of a more balanced use of the assigned resources. Employees are now used in the range 31%-77% of their time (whereas they were in the range 42%-48% with the previous assignment). The load for drones has slightly dropped, as a consequence of the better balanced load of employees and drones.

Search-Based Optimization

The analysis performed in the previous section requires the analysis of all possible combinations of numbers of resource instances. In the example in the previous section, the analysis was performed for the range [1..10] for both employees and drones. Considering more resources and wider ranges the analysis could be very costly.

It is a classical optimization problem and it can be solved by embedding the above-mentioned techniques for evaluating a specific resource assignment into some of the heuristic-based search optimization techniques. The possibility of using the techniques in the previous section into a simple gradient descent is illustrated next. Gradient descent is intended for finding local minimum (respectively, maximum) values starting from a specific point. After analyzing the function in the neighbouring points, the one with greatest value change is identified. Since the goal is to find a local minimum, the smallest value is chosen. If none of the neighbours is smaller, then the actual point is a local minimum.

By using the gradient descent algorithm on the running example, starting from the assignment (employee → 4, drone → 4), and using the optimization function and values in Section 4.3, the optimal solution can be found by following the sequence (4, 4) -→ (3, 5) -→ (2, 6) -→ (2, 7) -→ [START_REF] Arévalo | A metamodel to integrate business processes time perspective in BPMN 2.0[END_REF][START_REF] Christiansen | Formal semantics and implementation of BPMN 2.0 inclusive gateways[END_REF]. The values for which the optimization function is computed by our search function and the corresponding evaluations is depicted in Table 3. When compared to the analysis in Section 4.3, only 27 combinations were analyzed, instead of the 100 combinations analyzed above. More importantly, this dynamic search enables the exploration of greater search spaces. However, notice that this is a local minimum. A more sophisticated search algorithm may be required to find a global minimum, such as annealing or genetic algorithms.

Additional Examples

The BPMN process depicted in Figure 11 describes the recruitment process in an enterprise. This process shows how a candidate must fill in a hiring form, carry on a medical checkup, and in some cases apply for visa. Once the documentation is submitted, it is checked by the human resources office, which can decide to accept, reject or request additional documentation. If accepted, the candidate is informed, an assistant is in charge of preparing a welcome kit, and the technical staff is in charge of including the corresponding data in the enterprise's DB.

In this example, there are three resources modeled, namely human resources, the assistants, and the technical staff. The execution of the gradient descent algorithm on the recruitment example, starting from the assignment (hr → 3, assistant →

Related Work

Resource allocation-including optimization, allocation constraints, and optimal schedule allocation-has been extensively explored for the business process domain from several formal methods approaches. Schömig and Rau [START_REF] Schömig | A Petri Net Approach for the Performance Analysis of Business Processes[END_REF] use colored stochastic Petri nets to specify and analyze business processes in the presence of dynamic routing, simultaneous resource allocation, forking/joining of process-control threads, and priority-based queueing. In their work, each resource is equipped with properties grouped in a role defining if the resource is eligible to perform a certain activity. Li et al. [START_REF] Li | Performance Modeling and Analysis of Workflow[END_REF] introduce multidimensional workflow nets to model and analyze resource availability and workload. Oliveira et al. [START_REF] Oliveira | Quantitative Analysis of Resource-Constrained Business Processes[END_REF] use generalized stochastic Petri nets for correctness verification and performance evaluation of business processes. In their work, an activity can be associated to multiple roles and the completion of an activity can use a portion of the resources available to a role. They also propose metrics for evaluating process performance such as: the minimum number of resources needed for a role in order to complete Figure 12: BPMN process of the visa example a process, the expected number of activity instances when completing a process under the assumption of sufficient resources, and the expected activity response time.

The problem of understanding how bounded resources can impact the behavior of a process is approached by Netjes et al. by using colored Petri nets [START_REF] Netjes | Analysis of Resource-Constrained Processes with Colored Petri Nets[END_REF]. They introduce the notion of "flexible resource allocation" as a way to assign resources associated to a given role based on priorities. Havur et al. study the problem of resource allocation in business processes management systems where resources can be assigned constraints (e.g., time of availability) and have dependencies. Their technique is based on the answer set programming formalism and is capable of deriving optimal schedules. More recently, Sperl et al. [START_REF] Sperl | Resource utilization prediction in decision-intensive business processes[END_REF] describe a stochastic method for quantifying resource utilization relative to structural properties of the processes and historical executions.

The work presented in this paper can encompass most of the specification and analysis features in [START_REF] Schömig | A Petri Net Approach for the Performance Analysis of Business Processes[END_REF][START_REF] Li | Performance Modeling and Analysis of Workflow[END_REF][START_REF] Oliveira | Quantitative Analysis of Resource-Constrained Business Processes[END_REF][START_REF] Braghetto | Performance evaluation of business processes through a formal transformation to SAN[END_REF][START_REF] Sperl | Resource utilization prediction in decision-intensive business processes[END_REF] for BPMN. The rewriting logic specification presented here can be used to specify and verify resource availability and workload, several related measures, and the impact of bounded resources on a process behavior. Despite the fact that the main focus of this paper has been on resource allocation analysis and its optimization, the underlying rewriting logic semantics can simultaneously be a testbed for integrating and combining several types of mechanical formal analysis techniques such as algorithmic-based LTL model checking and statistical verification. Features such as multiple roles associated to an activity, flexible resource allocation, dependencies among the resources, and the use of historical data can be extensions of the present work.

There is also significant effort in the research community aimed at providing formal semantics and verification techniques for business processes using different formalisms, with Petri nets and process algebra as the preferred formalisms by most of them (see, e.g., [START_REF] Bruni | Applying process analysis to the italian egovernment enterprise architecture[END_REF][START_REF] Dijkman | Semantics and analysis of business process models in BPMN[END_REF][START_REF] Raedts | Transformation of BPMN models for behaviour analysis[END_REF]). [START_REF] Bruni | Applying process analysis to the italian egovernment enterprise architecture[END_REF] presents a business process monitoring platform that relies on the ProM mining tool and Petri nets tools for verifying synchronization, sojourn and waiting times. Raedts et al. [START_REF] Raedts | Transformation of BPMN models for behaviour analysis[END_REF] translate BPMN processes into Petri nets and use model checkers to analyze invariants such as deadlock freedom. The authors of [START_REF] Wong | A process semantics for BPMN[END_REF] present a formal semantics for BPMN by an encoding into the process algebra CSP. They show in [START_REF] Wong | Verifying business process compatibility[END_REF] how such a semantics can be used to verify compatibility between business participants in a collaboration. There is a further extension in [START_REF] Wong | A relative timed semantics for BPMN[END_REF] to propose a timed semantics of BPMN with delays. In [START_REF] Capel | Automating the transformation from BPMN models to CSP+T specifications[END_REF][START_REF] Morales | A formalization proposal of timed BPMN for compositional verification of business processes[END_REF][START_REF] Mendoza | Conceptual framework for business processes compositional verification[END_REF], the focus is on the semantics formalized in [START_REF] Wong | A process semantics for BPMN[END_REF][START_REF] Wong | A relative timed semantics for BPMN[END_REF] by proposing an automated transformation from BPMN to timed CSP, as well as composition verification techniques for checking properties using the FDR2 model checker. Poizat et al. [START_REF] Poizat | Checking business process evolution[END_REF] present an encoding of an untimed subset of BPMN into the LNT process algebra for supporting the analysis of process evolution. [START_REF] Durán | Verifying timed BPMN processes using maude[END_REF] presents an early attempt to represent BPMN processes enriched with time features in Maude. This work was extended in [START_REF] Durán | Stochastic analysis of BPMN with time in rewriting logic[END_REF] to provide automated verification of stochastic properties such as expected processing and synchronization times. Compared to these related works, the executable specification in rewriting logic presented in Section 3 can be seen as a semantic framework for BPMN, yet it is not the primary goal of this paper. The main difference with respect to those works is that the focus here is in the formal specification and verification of quantitative aspects of processes and their resources.

Several works propose extensions of BPMN with time constructs, see, e.g. [START_REF] Arévalo | A metamodel to integrate business processes time perspective in BPMN 2.0[END_REF][START_REF] Gagné | Time-bpmn[END_REF]. In [START_REF] Gagné | Time-bpmn[END_REF], the authors present Time-BPMN, an extension of BPMN to represent various constraints and dependencies that may arise when modelling business processes. The temporal constraints are used to control the beginning/end of an activity, whereas temporal dependencies involve two activities and indicate the relationship between their respective beginning/end. In [START_REF] Arévalo | A metamodel to integrate business processes time perspective in BPMN 2.0[END_REF], a metamodel-based approach to integrate temporal constraints and dependencies is introduced. The time aspects are specified using rules and OCL constraints capture the semantics of these rules. [START_REF] Gagné | Time-bpmn[END_REF] introduces time as a first class citizen to the BPMN standard whereas [START_REF] Arévalo | A metamodel to integrate business processes time perspective in BPMN 2.0[END_REF] proposes an MDE-based approach that enrich the existing BPMN control-flow graph model. In comparison, the goal in the present paper was not to extend the BPMN notation to take expressive modeling of time constraints into account, but to rely on usual time representation (duration of task and flow) and concentrate on the analysis of key timing properties in these models.

The expected processing time, the expected synchronization time for merge gateways, and other similar measures have been considered by some authors to evaluate deployed processes by, e.g., collecting timestamps from logs (see, e.g., [START_REF] Van Der Aalst | Process Mining -Data Science in Action[END_REF][START_REF] Van Der Aalst | Discovering workflow performance models from timed logs[END_REF][START_REF] Bruni | Applying process analysis to the italian egovernment enterprise architecture[END_REF]). Task duration has been specified using literal values in [START_REF] Braghetto | Performance evaluation of business processes through a formal transformation to SAN[END_REF][START_REF] Bruni | Applying process analysis to the italian egovernment enterprise architecture[END_REF][START_REF] Oliveira | Modeling and analyzing resource-constrained business processes[END_REF] and probabilities have been used to specify the likelihood of branching in [START_REF] Braghetto | Performance evaluation of business processes through a formal transformation to SAN[END_REF][START_REF] Oliveira | Modeling and analyzing resource-constrained business processes[END_REF]. Inclusive gateways are only supported in [START_REF] Braghetto | Performance evaluation of business processes through a formal transformation to SAN[END_REF][START_REF] Mateescu | Quantifying the parallelism in BPMN processes using model checking[END_REF][START_REF] Poizat | Checking business process evolution[END_REF] and loops only in [START_REF] El-Saber | BPMN formalization and verification using maude[END_REF][START_REF] Oliveira | Modeling and analyzing resource-constrained business processes[END_REF][START_REF] Poizat | Checking the realizability of BPMN 2.0 choreographies[END_REF]. The goal in most of these works is, however, quite different. Overall, Oliveira et al. [START_REF] Oliveira | Modeling and analyzing resource-constrained business processes[END_REF] possibly have the closest approach to the one presented in this paper: task durations and probabilities for exclusive branching can be specified, and loops and continuous time are considered. However, they cannot specify durations with stochastic expressions and do not perform analysis on shared resources.

Conclusion

This paper presented a solution to the problem of formal modeling and mechanical analysis of resource allocation for BPMN processes. The approach relies on the rewriting logic semantic framework and uses the Maude system (and language) to stochastically simulate multiple concurrent executions of a process instance that compete for the shared resources. The encoding of the BPMN syntax and execution semantics in rewriting logic supports an expressive subset of BPMN consisting of: activity and collaboration diagrams, several types of gateways, timed flows and tasks, probabilities for exclusive and inclusive split gateways, unbalanced workflows, and looping behavior. A main contribution of this work is to demonstrate how several properties of interest on resource allocation can be automatically analyzed, including resource charge over time and usage percentage for each resource replica. This paper also showed how the optimal allocation of resources can be computed by formulating a multi-objective optimization problem. The approach was validated on several examples, including real-world processes with large workloads.

As far as future work is concerned, a first perspective is to move from design time to runtime. The idea would be to dynamically adjust the amount of resources depending on the evolution of the workload. It is worth noting that this runtime solution would apply in some specific cases only when resources are available and could be used on-demand. A second perspective aims at using the results of the analysis proposed in this paper to address refactoring issues of a process in order to make it more efficient. As an example, the analysis may emphasize the unnecessary use of certain gateways resulting in bottlenecks in the execution. Ideally, the ultimate goal would be to have such a refactoring process as automated as possible for simplifying this optimization task.

Figure 1 :

 1 Figure 1: Supported BPMN syntax.

Figure 2 :

 2 Figure 2: Running example: parcel delivery by drones.

Figure 3 :

 3 Figure 3: Running example: representation in Maude of the parcel delivery process.

Figure 4 :

 4 Figure 4: Representation of the Simulation object.

  Simulation | tokens : (token(TId, NId, 0) Tks), ... Atts1 > 6 < CId : Counter | counter : N > 7 => < PId : Process | nodes : (start(NId, FId), Nodes), 8 flows : (flow(FId, SE), Flows), 9 Atts > 10 < SId : Simulation | tokens : insert(Tks, token(TId, FId, T')), ... Atts1 > 11 < CId : Counter | counter : N' > 12 if {T', N'} := eval(SE, N) .

Figure 5 :

 5 Figure 5: Start event processing.

Figure 6 :

 6 Figure 6: Task initiation rule.

Figure 7 :

 7 Figure 7: Event-based split gateway: one message available.

Figure 8 :

 8 Figure 8: Running example (400 instances): average execution time.

Figure 9 :

 9 Figure 9: Running example (400 instances): average usage percentage per employee (left) and drone (right).

  depicts normalized costs and average execution times. If the preference is to minimize delivery time, specified by, say, weights w cost = 0.4 and w aet = 0.6, the best solution would be the one marked in the figure, representing the combination of 3 employees and 8 drones.

Figure 10 :

 10 Figure 10: Cost vs. average execution time.

Figure 11 :

 11 Figure 11: BPMN process of the recruitment example.

Table 1 :

 1 Num. AET Var AST g8 AST ee Experimental results for the running example (2 employees, 3 drones).

				Total	Resources		Anal.
	inst.			time GTU e GTU 1 e UP e GTU d GTU 1 d UP d time
	100 107 190 70	57	327 316 158 48 852	284	87	6s
	200 160 37	81	108 599 506 253 42 1665 555	93	30s
	400 301 213 107	252 1202 1202 601 43 3493 1164 97 225s
	800 550 4	156	501 2367 1978 989 42 6953 2318 98 1748s
	1600 910 56 250	862 4263 4005 2003 47 12648 4216 99 11828s

Table 2 :

 2 Experimental results for the running example (3 employees, 8 drones).

Acknowledgments. The first author was supported in part by projects PGC2018-094905-B-I00 (Spanish MINECO/FEDER) and UMA18-FEDERJA-180 (FEDER Andalucía). The second author has been supported in part by CAPES, Colciencias, and INRIA via the STIC AmSud project "EPIC: EPistemic Interactive Concurrency" (Proc. No 88881.117603/2016-01), and Capital Semilla 2017 project "SCORES: Stochastic Concurrency in Rewrite-based Probabilistic Models" (Proj. No. 020100610).

0,591 0,591 0,567 0,687 2 0,560 0,492 0,467 0,425 0,428 0,446 3 0,514 0,621 0,508 0,480 0,464 0,420 0,528 4 0,622 0,685 0,653 0,527 0,461 0,508 0,466 5 0,626 0,731 0,603 6 7 8 9 10 The BPMN process depicted in Figure 12 describes a visa application process. The visa requester must provide a scanned copy of her passport and pay certain fees. The scanned passport is checked by an employee. If everything is all right, the visa is granted and delivered to the requester.

In this example, there are two resources modeled, namely the employee and the visa printer. Note that the task Deliver visa requires an instance of each of the resources. The execution of the gradient descent algorithm on the visa application example, starting from the assignment (employee → 1, printer → 5), using the optimization function and values in Section 4.3, with costs (per hour) (employee → 50, printer → 5), a local optimal solution is found by following the sequence


David Poulard
email: david.poulard@ifsttar.fr.

François Bermond

Sabine Compigne

Karine Bruyère

In-vivo analysis of thoracic mechanical response under belt loading: the role of Body Mass Index in thorax stiffness

Keywords: Biomechanics, sled test, BMI, stereovision, thoracic stiffness, volunteer

Thoracic injuries are a major cause of mortality in frontal collisions, especially for elderly and obese people. Car occupant individual characteristics like BMI are known to influence human vulnerability in crashes. In the present study, thoracic mechanical response of volunteers quantified by optical method was linked to individual characteristics. 13 relaxed volunteers of different anthropometries, genders and age were submitted to non-injurious sled tests (4 g, 8 km/h) with a sled buck representing the environment of a front passenger restrained by a 3point belt. A resulting shoulder belt force was computed using the external and internal shoulder belt loads and considering shoulder belt geometry. The mid sternal deflection was calculated as the distance variation between markers placed at mid-sternum and the 7 th vertebra spinous process of the subject. Force-deflection curves were constructed using resulting shoulder belt force and midsternal deflection. Average maximum chest compression was 7.9 ± 2.3 % and no significant difference was observed between overweight subjects (BMI≥25kg/m²) and normal subject (BMI<25kg/m²). The overweight subjects exhibited significantly greater resultant belt forces than normal subjects (715 ± 132 N vs. 527 ± 111 N, p < 0.05), higher effective stiffness (30.9 ± 10.6 N/mm vs. 19.6 ± 8.9 N/mm, p < 0.05) and lower dynamic stiffness (42.7 ± 8.71 N/mm vs. 61.7 ± 15.5 N/mm, p < 0.05).

Introduction

Although significant improvements have been achieved in mitigating road traffic fatalities, frontal impacts play a predominant role in the frequency of road traffic fatalities as they account for up to 46% of the mortality [START_REF] Klanner | Status report and future development of the Euro NCAP program[END_REF]. Statistics show that in France thoracic injuries were the first cause of death among car occupants in frontal collisions and the main source of serious casualties of elderly occupants [START_REF] Lafont | Vieillissement et gravité des atteintes lésionnelles des victimes d'accident de la circulation routière: Ageing and injury severity among road traffic accident victims[END_REF][START_REF] Ndiaye | Thorax vulnerability of female car users reported by the road accident field[END_REF]. Moreover, if we focus on moderate injuries sustained by the elderly population, rib and sternum fractures are mainly observed [START_REF] Ndiaye | Thorax vulnerability of female car users reported by the road accident field[END_REF].

First studies on Post Mortem Human Surrogates (PMHS) focused on generating thoracic force-deflection corridors using a hub-impact test condition [START_REF] Kroell | Impact tolerance and response of the human thorax. biomechanics of impact injury and injury tolerances of the thorax-shoulder complex[END_REF][START_REF] Lobdell | Impact response of the human thorax[END_REF]Nahum A.M., 1975;[START_REF] Lau | The viscous criterion -bases and applications of an injury severity index for soft tissues[END_REF]. These corridors were used in the development of frontal impact dummies. A statistical relationship was established between midsternal compression and thorax injury risk [START_REF] Kroell | Impact tolerance and response of the human thorax. biomechanics of impact injury and injury tolerances of the thorax-shoulder complex[END_REF][START_REF] Mertz | Hybrid III sternal deflection association with thoracic injury severities of occupants restrained with force limiting belts biomechanics of impact injury and injury tolerances of the thorax-shoulder complex[END_REF].

However, force-deflection corridors of the thorax were found to be strongly dependent on the load distribution [START_REF] Kent | Thoracic response to dynamic, non-impact loading from a hub, distributed belt, diagonal belt, and double diagonal belts[END_REF], and particularly on the anatomical area that supports the load [START_REF] Shaw | Quasistatic and dynamic thoracic loading tests: cadaveric torsos[END_REF]. Moreover, thoracic mechanical response depends strongly on the geometrical characteristics of the ribcage and on its biological material properties [START_REF] Kent | Structural and material changes in the aging thorax and their role in crash protection for older occupants[END_REF][START_REF] Gayzik | Quantification of age-related shape change of the human rib cage through geometric morphometrics[END_REF]. Authors agree that the deformation of the rib cage under an anteroposterior loading is divided into rib rotation relative to the vertebrae and sternum, and rib deflection [START_REF] Kent | Structural and material changes in the aging thorax and their role in crash protection for older occupants[END_REF]. More accurately, it is assumed that according to the rib orientation, the applied force causes rib deformations and rotations at various degrees [START_REF] Kent | Structural and material changes in the aging thorax and their role in crash protection for older occupants[END_REF]. Specifically, for an individual with a "vertical" rib cage, an applied anteroposterior load produces rotations of the rib joints and bone deformation whereas, for individual with a "horizontal" rib cage, the load force acts straight in the plane of the ribs and implies primarily bone deformation. A link was found between the initial rib slope and the amount of rotation and deformation [START_REF] Vezin | Structural characterization of human rib cage behavior under dynamic loading[END_REF]. Thus, it seems that geometrical characteristics of the ribcage may predispose ribs to fracturing.

In addition, it was well established that tolerance of the human thorax under dynamic loading decreases as age increases [START_REF] Ndiaye | Thorax vulnerability of female car users reported by the road accident field[END_REF][START_REF] Kroell | Impact tolerance and response of the human thorax. biomechanics of impact injury and injury tolerances of the thorax-shoulder complex[END_REF][START_REF] Zhou | Age effects on thoracic injury tolerance[END_REF][START_REF] Kent | Structural and material changes in the aging thorax and their role in crash protection for older occupants[END_REF]. Physiological phenomena of ageing could be described as the combination of mechanical and structural modifications associated with a decrease in chest deflection tolerance. The elastic modulus and the tensile strength of human bones decrease with age [START_REF] Zhou | Age effects on thoracic injury tolerance[END_REF]. Similarly, a decrease in the tensile strength of hyaline rib cartilage and cortical thickness in the rib sections was observed [START_REF] Takahashi | Age and sex related changes in the amount of cortex of normal human ribs[END_REF][START_REF] Yamada | Strength of Biological Materials[END_REF][START_REF] Kemper | The biomechanics of human ribs: material and structural properties from dynamic tension and bending tests[END_REF]. The age-related changes in thorax shape were also studied and authors assumed that the ribs became more horizontal with age and the slope of ribs in sagittal plane decreased with age [START_REF] Oskvig | Special problems in the elderly[END_REF][START_REF] Kent | Structural and material changes in the aging thorax and their role in crash protection for older occupants[END_REF][START_REF] Gayzik | Quantification of age-related shape change of the human rib cage through geometric morphometrics[END_REF]). Yet, some authors also found a correlation between BMI and the thorax shape [START_REF] Berthet | Assessment and analysis of the human rib lateral slopes[END_REF][START_REF] Gayzik | Quantification of age-related shape change of the human rib cage through geometric morphometrics[END_REF]. Body Mass Index (BMI) was associated with increased risk of mortality and increased risk of severe injury with a predominant occurrence for rib fractures and pulmonary contusions [START_REF] Boulanger | Body habitus as a predictor of injury pattern after blunt trauma[END_REF][START_REF] Mock | The relationship between body weight and risk of death and serious injury in motor vehicle crashes[END_REF]. Obese subjects were linked to unfavorable kinematics (greater head and pelvis excursions) and may be linked to higher chest compression [START_REF] Forman | The effect of obesity on the restraint of automobile occupants[END_REF]. In the vast majority of developed countries, the prevalence of overweight and obesity in older adult population is very high [START_REF] Gutierrez-Fisac | Prevalence of overweight and obesity in elderly people in Spain[END_REF].

Thus, in order to study the variability of the population in age, anthropometry and for both genders, in vivo experiments seem necessary. Except early in-vivo experimental studies [START_REF] Armstrong | Human muscular restraint during sled deceleration[END_REF][START_REF] Lobdell | Impact response of the human thorax[END_REF][START_REF] Stalnaker | Human torso response to blunt trauma[END_REF] which brought knowledge on the human body tolerance, most of them were used to analyze the physiology influence on the body response to dynamic loading -muscle tone [START_REF] Patrick | Impact force-deflection of the human thorax[END_REF][START_REF] Backaitis | Chest deflection characteristics of volunteers and hybrid III dummies[END_REF][START_REF] Kemper | Effects of pre-impact bracing on chest compression of human occupants in low-speed frontal sled tests[END_REF]-or to assess the response of specific populations such as children for non-injurious levels [START_REF] Sandoz | Mechanical properties assessment of child trunk[END_REF]Arbogast et al., 2009a;Arbogast et al., 2009b).

In the present paper, we present the analysis of thoracic mechanical response under belt loading of adults of various anthropometries and age subjected to a low deceleration pulse [START_REF] Poulard | In-vivo test protocol for analysis of thorax response under belt loading[END_REF]. In particular, the influence of overweight (BMI≥25kg/m²) on thoracic mechanical response will be assessed.

Materials and methods

This in vivo protocol was reviewed and approved by the French ethical committee Comité de Protection des Personnes Sud-Est II at Lyon in February 2011.

Test device

Ifsttar shock sled was set up with a standard sedan car environment of a front passenger restrained by a 3-point belt (Fig. 1). Posterior support consisted of a belt band tensed between the two backseat rods. Weight was added to the sled according to the weight of the subject (up to 110 kg) in order to maintain a constant sled mass (640 kg). The propelling system of the sled was chosen to allow low impact velocities with a good repeatability in pulses and to limit the maximum impact speed; it was made of two rubber bands of diameter 20 mm mounted on each side of the sled.

Deceleration pulse

As described in [START_REF] Poulard | In-vivo test protocol for analysis of thorax response under belt loading[END_REF], the deceleration pulse was chosen between 2 pulses measured on 2 public demonstrators and it was found similar to the one of a recent study (Arbogast et al., 2009a). The sled was pulled rearward by an electrical winch until reaching 165% of the two rubber bands' initial length so that an impact speed of 8 km/h was reached.

The deceleration system was obtained by laminating a polyurethane tube by a spear ended in an olive shape (hardness 95 SHA, 46.5 mm diameter spear). The sled was stopped with a maximum deceleration of 4 g during 120 ms.

Tested subjects

13 volunteers were tested. They were between 19 and 65 years of age and consisted of both males and females. They were chosen on criteria of height based upon anthropometry charts [START_REF] Jürgens | International data on anthropometry[END_REF]) and various BMI. Volunteer characteristics are listed in Table 1.

For physiological issues, subjects with existing endocrine disorders, mediastinal pathologies, previous or current diseases of the head, neck, spine were excluded from the tests. Only volunteers over 40 years of age with no osteopenia or osteoporosis (total femoral and lumbar Tscore>-1) were included in sled testing experiments. Prior to the testing, a medical doctor conducted an examination of each subject to confirm eligibility.

Two sled tests were performed on each volunteer. For each test, they were asked to adopt a relaxed posture.

Instrumentation

The longitudinal deceleration of the sled was recorded by a sled-mounted uniaxial accelerometer. Three belt webbing load cells were installed: on the external end (FB3) and on the internal end (FB4) of the shoulder belt and on the external side of the lap belt. Maximum value of FB3 and FB4 are noted FB3 max and FB4 max respectively. Reaction forces were measured by four uniaxial load cells placed under each footrest and under the seat pan. Data were acquired at 20 kHz during 5 s and filtered using CFC Butterworth filters according to the standard norm SAE J211.

Stereovision

Markers were placed on volunteers at anatomical landmarks:

-Spine: spinous processes of C7, T1, T7, and L4 and all along the spine, -Torso : on left and right acromion and clavicles, jugular notch (JN), mid-sternum (MS), xiphoid process (XP) and all along the sternum. Additional landmarks were additionally placed 5 cm and 10 cm to the left and right side of chest and on the inferior contour of the ribcage.

In addition, markers were placed on the sled in order to follow sled displacement.

Four high speed cameras were positioned around the sled for recording at 1000 fps the position of markers placed on subjects and sled (Fig. 2). Time zero was given by a flash visible by all cameras and synchronised with the data acquisition system. After the tests, the scene was calibrated by using a calibration object and the three dimensional coordinates of the subject markers were computed using the Direct Linear Transformation method also called stereovision [START_REF] Abdel-Aziz | Direct linear transformation from comparator coordinates into object space coordinates in close-range photogrammetry[END_REF]. Accuracy of the optical method was assessed in a previous study [START_REF] Poulard | In-vivo test protocol for analysis of thorax response under belt loading[END_REF]. 3D trajectories were filtered at CFC-60 Butterworth.

Anatomical landmarks were obstructed from view during some frames during the test, and thus, the 3D reconstruction by stereovision was not possible during these frames. Therefore the missing part of the trajectory was interpolated using other landmarks available at the same frame (see Appendices for more details).

Data analysis

The mid sternal deflection (D) was calculated as the distance variation between markers placed at mid-sternum and at 7th vertebra spinous process of the subject. The initial zero of deflection was defined as the midsternal thickness measured after reaching 5% of the maximum value of external shoulder belt load (see midsternal thickness in Table 1).

Midsternal compression C was defined as the ratio of the mid sternal deflection to the midsternal thickness. Maximum values of D and C are noted D max and C max respectively.

Resulting

shoulder belt force (FRes) is defined as

) cos( 2 2 2          FB4 FB3 FB4 FB3 FRes
and Fig. 3 as a combination of the external shoulder belt load (FB3), the internal shoulder belt load (FB4) and the angle made by the belt (α) calculated using the extreme markers of the belt [START_REF] Eickhoff | Optimization of seat belt bucke motion for reducing chest deflection, using rib eye sensors[END_REF]. A timeinterpolation of FB3 and FB4 was performed for combining angles values and loads values.

Maximum value of FRes is noted FRes max .

FRes is a 2D approximation of the 3D resulting belt force acting on the thorax. Nevertheless, the obtained FRes force vector is correctly aligned with the mid sternal deflection (D) vector (see Appendices for more details).

Force-deflection curves were built using FRes and mid sternal deflection D. Dynamic Stiffness (DS) was computed as the interpolated ratio of the resulting shoulder belt load force FRes over the midsternal deflection D before FRes max was reached. Effective Stiffness (ES)

was computed from the ratio of FRes over D when D max was reached [START_REF] Kent | Structural and material changes in the aging thorax and their role in crash protection for older occupants[END_REF][START_REF] Sandoz | Mechanical properties assessment of child trunk[END_REF]. While DS represents the thorax behaviour during the beginning of the load, ES represents its behaviour without taking account of this part.

The statistical analysis was carried out with R software for Windows ver. 2.14.2. Numerical variables were given as average ± standard deviation (SD). Groups were compared with permuted Student's t-test. Comparisons were considered as significant for p < 0.05.

Results

Definition of an individual thoracic mechanical response

Two sled tests were performed for each volunteer. The force-deflection responses obtained in these two tests were found significantly different. The results of the test for which the highest Dmax was measured, was only considered. Thoracic mechanical characteristics are synthetized in Table 2. The initial shoulder belt angle α is noted αini.

Seat belt loads

Separate curves from overweight (BMI≥25 kg/m², n=5) and normal subjects (BMI<25 kg/m², n=8) for the external shoulder belt load FB3 and the resulting shoulder load FRes are shown in Fig. 4. At the beginning of the shock, 5% of FB3max was reached at 55 ± 5 ms (n=13). At an average time of 118 ± 8 ms (n=13), FB3 and FRes reached their peak. No significant difference (p=0.18) were found in FB3 max between overweight subjects and normal subjects.

Average FRes max was 723 ± 131 N (n=5) for overweight subjects, significantly higher (p<0.02) than 535 ± 111 N (n=8) for normal subjects.

Midsternal compressions

Time histories of mid-sternal compressions for the overweight group (a) and the normal group (b) are shown in Fig. 5. At an average time of 136 ± 15 ms, the midsternal compression reached its peak, C max . There are no significant differences between overweight subjects (BMI≥25 kg/m², n=5) and the subjects in the normal group (BMI<25 kg/m², n=8) in C max (p=0.41).

The pulse applied in the present study is similar to the one applied by [START_REF] Kemper | Effects of pre-impact bracing on chest compression of human occupants in low-speed frontal sled tests[END_REF] on 5 males volunteers close to 50th percentile submitted to sled test. Average corridor of mid-sternal compressions are compared to those obtained in this study [START_REF] Kemper | Effects of pre-impact bracing on chest compression of human occupants in low-speed frontal sled tests[END_REF] on relaxed subjects (n=5) in Fig. 6. The average C max are similar in both studies (7.9% vs. 8%) even if midsternal compressions in the present study rose and dropped slower due to the difference in deceleration pulse characteristics, (4 g, 8 kph vs. 5 g, 9.7 kph).

Force-deflection characteristics

Fig. 7 shows the relationship between the resultant shoulder belt load FRes and the midsternal deflection D for the normal group (a) and the overweight group (b). The overweight subject thorax was more compliant at the beginning of the loading and FRes increases with deflection continuing on up to D max . There are no significant differences between overweight subjects (BMI≥25 kg/m², n=5) and the subjects in the normal group (BMI<25 kg/m², n=8) in D max (p=0.17) and as written before, FRes max was significantly higher for the overweight group than for normal group.

The average value of the dynamic stiffness (DS) was 54.8 ± 15.7 N/mm (n=13), greater than the average value of the effective stiffness (ES), 23.6 ± 10.1 N/mm (n=13). The different stiffness values are given in Table 2. Average ES was 30.9 ± 10.6 N/mm (n=5) for overweight subjects, significantly higher (p<0.02) than 19.6 ± 8.9 N/mm (n=5) for normal group, Fig. 8a.

Inversely, average DS was 42.7 ± 8.71 N/mm (n=5) for the overweight group, significantly lower (p<0.05) than 61.7 ± 15.5 N/mm (n=5) for the normal group, Fig. 8b. As only one elderly subject is available in the present study, no assumption can be made about the influence of age on thoracic stiffness.

Discussion

13 volunteers of different anthropometries, genders and age were submitted to a low deceleration pulse. Thoracic mechanical response quantified using non-invasive techniques was linked to individual characteristics and compared to literature. In particular, the influence of overweight on thoracic mechanical response was evaluated.

Thoracic mechanical responses obtained in the present protocol are similar to those obtained in sled tests [START_REF] Kemper | Effects of pre-impact bracing on chest compression of human occupants in low-speed frontal sled tests[END_REF]. The variability of thoracic mechanical response is higher in the present study than in [START_REF] Kemper | Effects of pre-impact bracing on chest compression of human occupants in low-speed frontal sled tests[END_REF] where only 50 th percentile males were considered.

The overweight volunteer group shows significantly higher ES and lower DS than the normal group. Increased compliance due to increased superficial tissue depth may cause the decrease in DS observed for overweight subjects at the beginning of the load. Since no significant difference was observed in D max in both group, the FRes max difference due to inertia, should explain consequently the higher ES observed for overweight subjects. Contrary to the trend highlighted by [START_REF] Forman | The effect of obesity on the restraint of automobile occupants[END_REF], no significant difference was found in C max for overweight group or for the few obese subjects (BMI≥30 kg/m², n=1) available in the present study.

As FRes max was also significantly linked to individual characteristics while FB3 max not; the implementation of loading direction in the computation of belt force resultant is important.

Although FRes is a 2D simplified computation, this approximation seems very adequate for the purposes of the current study and is easily reproducible for further studies.

Volunteer experiments overcome limitations related to selection criteria of human surrogates but have to face other challenges such as impact severity far below crash conditions and lower measurement accuracy due to the use of non-invasive techniques. Despite these limitations, this study could help in investigating the role of BMI in thoracic mechanical response variability and in assessing the thoracic mechanical response of a specific population, such the overweight and obese.

To the authors' knowledge, this is the first study to compare thoracic mechanical responses from overweight subjects to normal subjects using frontal sled tests with volunteers. The overweight subjects exhibited significantly greater resultant belt forces, higher effective stiffness and lower dynamic stiffness. Increased soft tissue depth and increased thorax may explain these observed differences.

The recruitment of volunteers has been continued including the elderly. It should allow isolating the influence of age on thoracic mechanical response independently from anthropometrical variations.

Fig. 1 .Fig. 2 .

 12 Figure legends:

Fig. 3 .

 3 Fig.3. Computation of the resulting shoulder belt load FRes as defined by[START_REF] Eickhoff | Optimization of seat belt bucke motion for reducing chest deflection, using rib eye sensors[END_REF].

Fig. 4 .

 4 Fig. 4. Time history of the seat belt loads, external shoulder belt force F B3 and the resultant shoulder belt force F Res for normal group (a,c) and overweight group (b,d).

Fig. 5 .

 5 Fig. 5. Time history of midsternal compressions for normal group (a) and overweight group (b).

Fig. 6 .

 6 Fig. 6. Comparison of midsternal compression corridors from the present study and Kemper et al. (2011).

Fig. 7 .

 7 Fig. 7. Force-deflection of thorax volunteers for normal group (a) and overweight group (b).

Fig. 8 .

 8 Fig. 8. Box plots of the effective stiffness (a) and dynamic stiffness (b) in N/mm for normal subjects and overweight subjects.

Fig. 9 .

 9 Fig. 9. Comparison of different methods of interpolations of mid-sternum (MS) landmark trajectory for the computation of mid-sternum deflection.

Fig. 10 .

 10 Fig. 10. Time history of the resulting shoulder belt load FRes projected towards the direction of the mid-sternal deflection vector (a) and its contribution to FRes magnitude (b).

Fig

  Fig. 9.

Table 1 .

 1 Anthropometric characteristics of volunteers.

Table 2 .

 2 Thoracic mechanical characteristics obtained on volunteers.

	Table 1
	ID	Gender	Age (yo)	Height (cm)	Weight (kg)	BMI (kg/m²)	Midsternal thickness (mm)
				Normal subjects
	200	F	22	176	58	19	179
	607	M	32	164	57	21	223
	862	F	41	161	58	22	224
	806	M	21	196	88	23	217
	869	M	19	167	64	23	226
	733	M	34	198	92	23	233
	957	F	23	155	58	24	187
	225	M	32	185	83	24	210
	Average	-	28	170	70	22	212
	SD	-	8	16	15	2	19
				Overweight subjects
	329	M	38	173	74	25	213
	948	M	25	169	73	26	225
	118	M	65	167	74	27	234
	442	M	35	176	86	28	258
	370	F	20	173	90	30	219
	Average	-	37	172	79	27	230
	SD	-	20	4	9	2	17
				All subjects
	Average	-	31	174	73	24	219
	SD	-	13	13	13	3	20
	Min	-	19	155	57	19	179
	Max	-	65	198	92	30	258

Acknowledgements

The authors would like to acknowledge all the volunteers who participated in this study. The authors would also like to acknowledge Patrick Joffrin, Alain Maupas, Jean-Luc Russo, Sophie Serindat, Loretta Leoni-Duplessy, Dr Georges Baraton for their contribution in the sled tests preparations and the Desgenettes Hospital for providing bone density examinations.

In addition, the authors would like to acknowledge Mitsutoshi Masuda from Toyota Motor Corporation for his advice.

Conflict of interest statement

This study was partly funded by Toyota Motor Europe. S. Compigne is employee of Toyota Motor Europe.

Appendices

Assessment of methods for the interpolations of mid-sternum landmark trajectory

Two methods were investigated to estimate the missing part of a marker trajectory. First method was a time-based interpolation, based upon cubic interpolation of the anterior and posterior frames available. Second method was a geometrical-based interpolation, based on dual-kriging method [START_REF] Krige | A statistical approach to some basic mine valuation problems on the Witwatersrand[END_REF] using other landmarks available at the same frame.

The two methods were assessed on a known mid-sternum (MS) trajectory during which 30 successive frames were deliberately hidden. The reconstructed MS trajectories were directly used in the computation of mid-sternal deflection. Fig. 9 shows the comparison of the sternal deflections computed from the reference MS trajectory and from the interpolated ones. The second method appeared more suitable since numerous successive frames interpolation could be applied without loss of the deflection peak pattern.

Assessment of FRes as a 2D approximation of the 3D belt force acting on the thorax

FRes provides the resultant force from two known forces with known angle between in a plane. It is a simplified computation of 3D resulting belt force acting on the thorax.

We investigated how the FRes force vector was aligned with the mid sternal deflection (D) vector used to construct the force-deflection curve. Time history of the resulting shoulder belt load FRes projected towards the direction of the mid sternal deflection vector (a) and its contribution to FRes magnitude (b) are shown in Fig. 10. Average contribution of the mid sternal deflection direction to FRes magnitude during the loading (dotted lines) was 94.2 ± 2.1 %. The 2D approximation used for FRes appeared very adequate for the purposes of the current study since the obtained FRes force vector was correctly aligned with the mid sternal deflection vector.