Julien Rey

Inès Rasolohery

Pierre Tufféry

Frédéric Guyon
email: frederic.guyon@univ-paris-diderot.fr

Gautier Moroy
email: gautier.moroy@univ-paris-diderot.fr

Pierre Tuff Éry

Fréd Éric Guyon

PatchSearch: a web server for off-target protein identification

published or not. The documents may come

INTRODUCTION

During the drug discovery process, binding sites comparison can assist in the identification of interactions of drugs with undesired targets (off-targets) and the understanding of adverse effects. Binding site comparison is also helpful for drug repositioning and ligand selectivity optimization. Consequently, different approaches have been developed for this purpose and include ligand-based and structure-based approaches [START_REF] Vulpetti | Chemogenomics in drug discovery: computational methods based on the comparison of binding sites[END_REF][START_REF] Jalencas | Identification of similar binding sites to detect distant polypharmacology[END_REF].

When based on the knowledge of the structure, off-target binding site identification faces the issue of structural plasticity, which hampers the identification of undesired binding partners. Different strategies have been considered and are mostly based on the fact that similar structures or regions of structure accessible to the solvent can be expected to bind similar ligands. Alignment-free methods perform an overall comparison of global properties and characteristics of binding sites such as shape, surface descriptors and physicochemical residue properties combined with atom types [START_REF] Ritchie | Fast computation, rotation and comparison of low resolution spherical harmonic molecular surfaces[END_REF][START_REF] Morris | Real spherical harmonic expansion coefficients as 3D shape descriptors for protein binding pocket and ligand comparisons[END_REF][START_REF] Hoffmann | A new protein binding pocket similarity measure based on comparison of clouds of atoms in 3D: application to ligand prediction[END_REF], Patch-Surfer [START_REF] Sael | Binding ligand prediction for proteins using partial matching of local surface patches[END_REF][START_REF] Sael | Detecting local ligand-binding site similarity in nonhomologous proteins by surface patch comparison[END_REF], PocketMatch [START_REF] Yeturu | PocketMatch: a new algorithm to compare binding sites in protein structures[END_REF], PocketFeature [START_REF] Liu | Using multiple microenvironments to find similar ligand-binding sites: application to kinase inhibitor binding[END_REF]. On the other hand, sequence order-independent alignments of residues or atoms are in general far more difficult to compute than alignment-free comparisons, but these methods allow for the identification of atoms or residues involved in the binding with a ligand. These methods are based on geometric hashing: TESS [START_REF] Wallace | TESS: a geometric hashing algorithm for deriving 3D coordinate templates for searching structural databases. Application to enzyme active sites[END_REF], SitesBase [START_REF] Gold | Fold independent structural comparisons of protein-ligand binding sites for exploring functional relationships[END_REF], SiteEngine [START_REF] Shulman-Peleg | Protein-protein interfaces: recognition of similar spatial and chemical organizations[END_REF] and I2I-SiteEngine [START_REF] Shulman-Peleg | Recognition of functional sites in protein structures[END_REF], MultiBind [START_REF] Shatsky | The multiple common point set problem and its application to molecule binding pattern detection[END_REF][START_REF] Shulman-Peleg | MultiBind and MAPPIS: webservers for multiple alignment of protein 3D-binding sites and their interactions[END_REF] and PCalign [START_REF] Cheng | PCalign: a method to quantify physicochemical similarity of protein-protein interfaces[END_REF], or on the Hungarian algorithm eMatchSite [START_REF] Brylinski | eMatchSite: sequence order-independent structure alignments of ligand binding pockets in protein models[END_REF]. A new approach based on deep learning has been recently published to compare binding site [START_REF] Pu | DeepDrug3D: classification of ligand-binding pockets in proteins with a convolutional neural network[END_REF]. Many methods also compute sequence order-independent alignment by searching for cliques in product graphs [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF]. The Bron-Kerbosh algorithm is the most efficient algorithm to search for all maximal cliques [START_REF] Bron | Algorithm 457: finding all cliques of an undirected graph[END_REF]. For this purpose, it is widely used, in particular in computational chemistry [START_REF] Gardiner | Graph-theoretic techniques for macromolecular docking[END_REF] and is recognized as being one of the most efficient in practice [START_REF] Cazals | A note on the problem of reporting maximal cliques[END_REF]. Many improved variants have since been described and more efficient algorithms for finding a maximum clique ex-W366 Nucleic Acids Research, 2019, Vol. 47, Web Server issue ists [START_REF] Cazals | A note on the problem of reporting maximal cliques[END_REF][START_REF] Konc | An improved branch and bound algorithm for the maximum clique problem[END_REF]. However, the Bron-Kerbosh algorithm provides a mean to explore all maximal cliques and therefore all possible matchings. The first methods developing this strategy have been applied to protein structure comparisons since early 90's [START_REF] Grindley | Identification of tertiary structure resemblance in proteins using a maximal common subgraph isomorphism algorithm[END_REF][START_REF] Gardiner | Clique-detection algorithms for matching three-dimensional molecular structures[END_REF], and more recently, clique algorithms have been used in CavBase [START_REF] Schmitt | A new method to detect related function among proteins independent of sequence and fold homology[END_REF] and eF-site [START_REF] Kinoshita | Identification of protein functions from a molecular surface database, eF-site[END_REF], SuMo [START_REF] Jambon | A new bioinformatic approach to detect common 3D sites in protein structures[END_REF]. PocketMatch, SiteEngine, eF-site, MultiBind and ProBis [START_REF] Konc | ProBiS algorithm for detection of structurally similar protein binding sites by local structural alignment[END_REF] are available as web servers (Supplementary Table S1). Most of the above approaches compare or align binding sites only. ProBis web server is the only one able to search for a binding site on the entire surface of proteins based on local structural alignments. ProBis web server requires a query structure of a protein-ligand complex. The user can select a query binding site which is compared to entries in the non-redundant PDB (nr-PDB) or to a user-supplied list of PDB identifiers.

Molecular docking approaches can be also used to identify protein target of a ligand and consequently help the detection of off-target protein. Thus, IdTarget web server was developed to predict possible binding targets of a small chemical molecule via a divide-and conquer docking approach [START_REF] Wang | idTarget: a web server for identifying protein targets of small chemical molecules with robust scoring functions and a divide-and-conquer docking approach[END_REF]. It requires an input ligand file for the target screening. The user can choose to perform the search of potential binding targets among two predefined datasets of PDB identifiers or a user-supplied list of PDB identifiers.

Recently, PatchSearch [START_REF] Rasolohery | PatchSearch: a fast computational method for off-target detection[END_REF] was developed to search for structurally conserved binding sites on the entire surface of a protein in order to help for the detection of potential off-target protein. It uses a quasi-clique approach which avoids a too stringent distance conservation between atoms and hence takes into account flexibility of binding sites. A quasi-clique is a dense subgraph. Our approach is similar to those used for dense subgraph or community detection in graph clustering [START_REF] Girvan | Community structure in social and biological networks[END_REF][START_REF] Clauset | Finding local community structure in networks[END_REF][START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. Cliques in correspondence graph involves the conservation of all internal distances between protein and patch surfaces. Based on Euclidean distance matrix properties, a well-chosen set of conserved internal distances is sufficient to ensure that all distances are equals or almost equals. Hence the complete node connection condition can be relaxed. In PatchSearch, a quasiclique is a correspondence subgraph including a clique with added nodes ensuring a spatial similarity. It is computed with a greedy algorithm which starts from cliques and at each step add nodes with at least four connections with the current quasi-clique and maximizing a 3D similarity score. The quasi-clique approach compared to a classical clique technique allows for fast detection of larger and new patches on protein surfaces which could potentially provide new yet unpredicted off-target binding sites.

PatchSearch has been validated against a number of widely used datasets [START_REF] Kahraman | Shape variation in protein binding pockets and their ligands[END_REF][START_REF] Hoffmann | A new protein binding pocket similarity measure based on comparison of clouds of atoms in 3D: application to ligand prediction[END_REF][START_REF] Gunasekaran | How different are structurally flexible and rigid binding sites? Sequence and structural features discriminating proteins that do and do not undergo conformational change upon ligand binding[END_REF][START_REF] Chen | Creation of a free, Internet-accessible database: the multiple target ligand database[END_REF][START_REF] Milletti | Predicting polypharmacology by binding site similarity: from kinases to the protein universe[END_REF][START_REF] Schomburg | Benchmark data sets for structure-based computational target prediction[END_REF] (Supplementary Table S2) with a large diversity of ligands. First, PatchSearch ability to recognize proteins binding the same ligand has been assessed on four datasets (Kahraman, Homogeneous, Gunasekaran and Milletti). Second, it was applied to three drugs (ursodeoxycholic acid, prazosin and naproxen) from Drugs/sc-PDB dataset and to three drugs (sunitinib, imatinib and sorafenib) used in cancer therapy from the Multiple Target Ligand Database. The aim of these experiments was to prove that PatchSearch manages to identify true offtarget proteins.

The main motivation of the present development of the PatchSearch web server is to provide the user with an automated workflow (Figure 1) to identify among a selection of proteins, which of them have a region sharing structural similarities with a ligand binding site. Through an energy minimization by steepest descent algorithm, the ligand conformation is locally optimized and scored within the binding site by smina program (41) using Vinardo scoring function [START_REF] Quiroga | Vinardo: A scoring function based on autodock vina improves scoring, docking, and virtual screening[END_REF].

MATERIALS AND METHODS

Patches and surfaces

A patch is a small piece of protein surface given by a set of solvent accessible atoms. Here, a patch is constituted by the solvent accessible atoms that have a distance smaller than 5 Å to the ligand. Each atom of the patch is assigned a label corresponding to its type: N, O, S, C and Ca (for carbons ␣). Aromatic rings are replaced by one or two pseudoatoms corresponding to one centroid for phenylalanine or tyrosine and two centroids for tryptophan). In this study, we are only interested in patch representing binding sites with ligands. The solvent accessibility is calculated using the program NACCESS (Hubbard, S.J. & Thornton, J.M. (1993), http://www.bioinf.manchester.ac.uk/naccess/) and all atoms with relative accessibility over 1% are conserved. This very small threshold allows to extract the largest surface involved in the interaction with the ligand.

For residues having at least one atom exposed, solvent accessible atoms and C␣ coordinates are stored for the structural similarities calculations.

Two patches are considered to be similar if it exists a matching associating equivalent atoms. The largest match is computed such that all inner distances between patch atoms are preserved in target surfaces within a given tolerance. This is obtained by computing a sequence orderindependent alignment, searching for cliques in product graphs [START_REF] Ullmann | An algorithm for subgraph isomorphism[END_REF] using the Bron-Kerbosh algorithm [START_REF] Bron | Algorithm 457: finding all cliques of an undirected graph[END_REF]. Yet, this clique strategy, though accurate, presents some drawbacks: comparisons of proteins with several thousand atoms can lead to some very large product graphs, and hence to a large amount of running time. In addition, a too stringent distance conservation between atom does not take into account flexibility of binding sites. Therefore, in Patch-Search a quasi-clique strategy is employed: a quasi-clique is constructed from a small core clique found with the Bron-Kerbosh algorithm with a distance tolerance of 1.2 Å, which is enriched using a greedy algorithm with less stringent distance tolerance set by default to 3 Å. The greedy algorithm optimizes the relative Binet-Cauchy (BC) structural scores between patch and matched atoms by merging cliques and adding atom matchings.

The BC score is a geometric correlation score [START_REF] Guyon | Fast protein fragment similarity scoring using a Binet-Cauchy kernel[END_REF]. The relative BC (rBC) score, is the BC score weighted by the percentage of retrieved atoms (N match ) relatively to the number of atoms of the patch (N patch ): r BC = N match N patch × BC

Rescoring step

For each targeted protein, the probable binding sites identified by PatchSearch are ranked according to the rBC score. [START_REF] Quiroga | Vinardo: A scoring function based on autodock vina improves scoring, docking, and virtual screening[END_REF]. During this step, the protein is rigid, the ligand is fully flexible and a maximum of 100 iterations of steepest descent algorithm are required.

The lowest resulting energy among the five minimized mod-els indicates which binding site can be considered as the most probable binding site at the protein surface.

WEB SERVER

The PatchSearch web server allows a local non-sequential searching for similar regions, called patches, on the entire protein surfaces, without any knowledge on binding site localizations or preliminary binding site detections. The service is fully embedded in the Mobyle framework [START_REF] Néron | Mobyle: a new full web bioinformatics framework[END_REF]. It embeds simple yet powerful data management features that allow the user to reproduce analyses. It gives to the users the possibility to create registered accounts, which allows user data and jobs to be maintained and managed across multiple work sessions and therefore to reuse and share data and results.

Input

There are two main inputs. First, the user must provide a PDB identifier containing at least one protein-ligand complex. The ligand, for which similar binding sites have to be detected, can be selected in the list of non-protein molecule in the PDB file. The selected ligand, the patch around the ligand and the protein surface can be visualized by NGL viewer [START_REF] Rose | NGL Viewer: a web application for molecular visualization[END_REF]. If the user knows the residue number and the chain of the ligand of interest, he can directly fill the fields in the form.

Second, the user must define a collection of structures in which the search will be performed. This can be done either on the form of a collection of proteins defined as a list of PDB identifiers provided by the user. At present, this list cannot exceed 5000 proteins. Alternatively, five predefined lists are available. They correspond to subsets of the PDB with a minimal sequence length of 50 residues. These subsets are created using the Advanced Search Interface in the Protein Data Bank server. Three lists correspond to human proteins at different sequence identities (30%: 7922 structures; 50%: 9118 structures; 70%: 9836 structures), a list of 14 191 eukaryotic proteins with 30% of sequence identity and a list of 14 650 prokaryotic proteins with 30% of sequence identity.

Output

The server returns patches similar to the query identified among the collection of proteins. Pairings between the query atoms and the atoms in the targeted protein surface are detailed and patches are scored by the relative Binet-Cauchy score. Similar patches are output to a table with the PDB identifier, the number of patch atoms, alignment length (the number of retrieved atoms), the RMSD between the query patch and the retrieved patch, the relative Binet-Cauchy score and docking values resulting from the Smina rescoring calculations.

All the alignment results, as well the estimated affinity and the binding mode of the ligand for each targeted protein are downloadable in separated files.

Patch visualization

An interactive page allows to browse the retrieved patches. The residues forming the retrieved patch are in 'lines' rep-resentation with C atoms in green, the ligand is in 'sticks' representation with C atoms in cyan and the targeted protein in white 'cartoon' representation. The best solution for each targeted proteins are ranked according to the affinity value. For the visualization, we use NGL viewer [START_REF] Rose | NGL Viewer: a web application for molecular visualization[END_REF] which takes advantage of WebGL capability of modern browsers for molecular graphics (Figure 2).

Execution times

Depending of the number of patch atoms and the protein size, average run times of searching one patch against an entire protein surface are <0.5 s (Supplementary Table S3). For the bigger predefined list of PDB containing 14 650 protein structures, the typical run times for a large patch formed by 100 atoms are of 1.5 h, but this depends on server load.

CASE STUDY

Distant off-targets

To illustrate PatchSearch effectiveness, we present an example of the identification of distant off-targets from the TOUGH-M1 dataset [START_REF] Govindaraj | Comparative assessment of strategies to identify similar ligand-binding pockets in proteins[END_REF]. The TOUGH-M1 dataset contains proteins with dissimilar global sequences and structures. The structural dissimilarity was measured by the TMscore [START_REF] Zhang | Scoring function for automated assessment of protein structure template quality[END_REF]. In this dataset, off-targets share no significant sequence (sequence identity < 30.0%) or fold similarity (TMscore < 0.4). This dataset is divided in two lists: a 'positive' list with proteins able to bind a chemically similar ligands and a 'negative' list in which the proteins interact with dissimilar ligands. Considering as input the binding site around Adenosine DiPhosphate (ADP) molecule in a myosin structure (PDB ID: 1lkx), we used the server to rank 20 proteins in 'positive' list and 20 proteins in 'negative' list for ADP. The patch around ADP was extracted from the 1lkx structure and compared to the entire surface of these 40 proteins. Note that all positive and negative proteins have very low sequence identities with the query --<23.0% and TM-scores <0.4. The results are reported in table 1 (Table 1).

The positive proteins have good affinity scores, with values less than -10.0 k cal/mol , indicating that a site on the surface has a good affinity for ADP. Therefore, the positive proteins are correctly detected as off-targets. For all negative proteins, except for 4nym, the binding affinity scores are poorer around -6 k cal/mol .

Results obtained for 4nym structure are interesting, because this protein is supposed to be unable to bind ADP. 4nym corresponds to the structure of GTPase HRas protein, which is involved in the activation of the Ras protein signal transduction. GTPase HRas protein is able to bind to Guanosine TriPhosphate (GTP) or Guanosine DiPhosphate (GDP). In this structure, GT-Pase HRas protein has been co-crystallized with Phosphoaminophosphonic Acid Guanylate Ester (ligand identifier for the PDB: GNP), a non-hydrolyzable analog of GTP, and the N-[1-(1H-indol-3-ylmethyl)piperidin-4-yl]-Ltryptophanamide (ligand identifier for the PDB: RND), a small molecule, altering experimentally the GTPase HRas activity. Logically, 4nym structure is considered as a negative protein for the binding site comparison between ADP The patch was extracted around ADP in the 1lkx myosin structure. Structural similarities between the patch and the entire surface of 40 PDB structures known to be able to interact, i.e. 'positive', or not, 'negative', with ADP or similar ligands, according to the TOUGH1-dataset. The sequence identity and the structure similarity were calculated between 1lkx against the each 40 PDB structures. The estimated affinity reported is the best score, computed by Smina program with Vinardo scoring function, between ADP and the five potential sites with the highest structural similarities according to the rBC score.

and RND, because these molecules are chemically dissimilar. However, the PatchSearch web server searches structural similarities onto the entire surface of the protein. The favorable computed binding affinity is basically due to the fact that PatchSearch has retrieved the GTP/GDP binding site of GTPase HRas protein (Figure 3). This unexpected result can be explained by the high chemical similarities between ADP and GDP. It reinforces strengthens the interest of taking into account the entire surface protein, not only comparing predetermined binding sites for the finding of new off-targets.

DISCUSSION AND CONCLUSION

Based on a new algorithmic approach, PatchSearch web server allows fast structural comparisons between a binding site and the entire protein surfaces of a user-supplied collection of protein structures. PatchSearch recognizes structural similarities based on geometry and physicochemical conditions by the matching of equivalent atoms. PatchSearch is able to perform the search on the entire protein surface to identify structural similarities without a priori knowledge about the binding sites. The use of quasi-cliques approach allows the detection of structurally flexible binding sites and the detection of similar binding sites having local structural distortions. In addition, PatchSearch web server is fast, the structural similarity calculations on ∼15 000 PDB structures are done in <2.0 h. To this respect, PatchSearch improves clearly over other methods dedicated only to the binding sites comparisons.

The PatchSearch web server benefits from a user-friendly submission and visualization interface. The transposed ligand within the potential binding site can be downloaded in PDB format for offline in-depth analysis. Consequently, the use of PatchSearch web server may be a preliminary step for the discovering of new interactions and hence is a valuable tool for predicting adverse effects, for helping in drug repositioning studies or for modifying a drug in a way that maintains binding to the intended target, but reduces binding to undesired proteins.

PatchSearch will often return a surprisingly high number of hits, in most cases with unknown status. However, it is important to keep in mind that the unwanted ligand interaction on a protein surface might show no effect on the protein behavior, especially if the potential binding site is not implicated in the protein biological activity, like enzymatic reactions or the interactions with partners. The potential off-targets have to be ascertained by further computational or experimental analysis.

Future directions are to extend PatchSearch to larger binding sites, such as protein-protein or protein-peptide binding sites, that have gain in recent years increasing interests.

Figure 1 .

 1 Figure 1. PatchSearch web server flowchart.

Figure 2 .

 2 Figure 2. Visualization of hits in the PatchSearch web server. (A) results table of proteins with patches similar to a query patch. (B) NGL viewer displaying the retrieved patch for 1efp structure. The protein is in white cartoon, ligand in thick sticks with C atoms in cyan and the residues of the retrieved patch are represented in lines with C atoms in green.

Figure 3 .

 3 Figure 3. The position of ADP, in sticks with C atoms in magenta, proposed by the PatchSearch web server and the experimental position of GNP, in sticks with C atoms in green, in the 4nym structure.

Table 1 .

 1 Example of PatchSearch results: identification of distant off-targets for ADP

			Structure
		Sequence identity	similarity	TOUGH1-dataset	Estimated Affinity
	PDB ID	(%)	(TM-score)	list	(k cal/mol )
	1fqj	20.5	0.398	positive	-12.934
	1in6	20.0	0.333	positive	-11.977
	2o8b	22.1	0.270	positive	-11.923
	1shz	18.5	0.386	positive	-11.853
	4zkd	21.3	0.298	positive	-11.537
	2fna	22.0	0.341	positive	-11.510
	4nym	18.1	0.282	negative	-11.390
	4d25	21.2	0.327	positive	-11.387
	3fwy	23.5	0.365	positive	-11.361
	1d2e	21.1	0.284	positive	-11.130
	5bn3	21.6	0.271	positive	-10.829
	1tq6	19.9	0.329	positive	-10.809
	3iev	23.2	0.352	positive	-10.692
	3u5z	21.3	0.299	positive	-10.616
	1kk3	20.4	0.297	positive	-10.532
	4kxf	19.3	0.298	positive	-10.491
	1dg1	16.3	0.295	positive	-10.444
	4djt	22.3	0.376	positive	-10.378
	1sxj	22.9	0.307	positive	-10.368
	3r7w	20.6	0.326	positive	-10.290
	4yj1	15.1	0.279	positive	-10.271
	2ozr	13.4	0.346	negative	-6.907
	3d3h	22.3	0.355	negative	-6.861
	1tve	18.8	0.339	negative	-6.780
	5ai9	20.4	0.278	negative	-6.572
	2vax	19.6	0.340	negative	-6.484
	4rzm	18.1	0.339	negative	-6.425
	1dsy	21.7	0.372	negative	-6.416
	5hes	17.9	0.322	negative	-6.349
	1y0g	19.0	0.343	negative	-6.328
	2rjp	16.9	0.335	negative	-6.319
	2f9a	20.6	0.294	negative	-6.298
	3n0t	19.3	0.348	negative	-6.223
	1kr1	18.4	0.370	negative	-6.214
	1fbo	19.4	0.297	negative	-6.152
	3njj	22.1	0.371	negative	-6.139
	3od2	19.3	0.362	negative	-6.123
	4s3r	20.8	0.283	negative	-6.059
	2zjf	20.5	0.350	negative	-6.039
	1lqy	19.6	0.361	negative	-5.923

FUNDING

Agence Nationale de la Recherche ANR-IA-2011-IFB [ANR-11-INSB-0013]. Funding for open access charge: Université Paris Diderot; INSERM [U1133].

DATA AVAILABILITY

The server is freely available via a user-friendly web interface at: http://bioserv.rpbs.univ-paris-diderot.fr/services/ PatchSearch/

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online. Conflict of interest statement. None declared.


Yoro Diouf
email: yorodiouf@gmail.com

Mamadou Oumar Maiga

Kaba Traore

Oumar Y Maïga

Mamadou K Traore

A Theoretical approach to the computational complexity measure of abstract DEVS simulators

Keywords: Discrete-Event Systems Specification (DEVS), Modeling and Simulation (M&S), computational complexity, simulation algorithm. Mathematics Subject Classification 2020: 00A72

DEVS is a sound Modeling and Simulation (M&S) framework that describes a model in a modular and hierarchical way. It comes along with an abstract simulation algorithm which defines its operational semantics. Many variants of such an algorithm have been proposed by DEVS researchers. Yet, the proper interpretation and analysis of the com-putational complexity of such approaches have not been systematically addressed and defined. As systems become larger and more complex, the efficiency of the DEVS simu-lation algorithms in terms of time complexity measure becomes a major issue. Therefore, it is necessary to devise a method for computing this complexity. This paper proposes a generic method to address such an issue, taking advantage of the recursion embedded in the triggered-by-message principle of the DEVS simulation protocol. The applicability of the method is shown through the complexity analysis of various DEVS simulation algorithms.

Introduction

DEVS (Discrete-Event Systems Specification) is widely recognized as a universal framework for discrete-event simulation. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] It proposes a sound formalism that comes along with an abstract simulator, which defines its execution semantics. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] An abstract DEVS simulator is a technology-agnostic description of a correct DEVS simulation engine. Therefore, it abstracts away all implementation details and only focuses on the description of the operations to be performed in order to ensure the correct generation of DEVS models' behavior. This abstract description is based on a hierarchical architecture of abstract components called simulators and coordinators, and a generic message-based inter-component communication scheme called the DEVS simulation protocol. In this paper, the terms "DEVS abstract simulator" and "DEVS simulation algorithm" are interchangeably used. Since the introduction of the original abstract DEVS simulator, 1 several variants have been proposed in the literature, including sequential simulation algorithms [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF][START_REF] Chow | Parallel DEVS: a parallel, hierarchical, modular modelling formalism and its distributed[END_REF][START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF][START_REF] Schwatinski | An advanced simulation approach for parallel DEVS with ports[END_REF][START_REF] Himmelspach | Sequential processing of PDEVS models[END_REF] as well as parallel and distributed versions. [START_REF] Wainer | Discrete Event Modelling and Simulation: A Practitioner's Approach[END_REF][START_REF] Syriani | Modelling and simulation-based design of a distributed DEVS simulator[END_REF][START_REF] Bisgambiglia | DecPDEVS: New simulation algorithms to improve message handling in PDEVS[END_REF][START_REF] Martin | A new simulation algorithm for PDEVS models with time advance zero[END_REF][START_REF] Cárdenas | A DEVS simulation algorithm based on shared memory for enhancing performance[END_REF] Yet, the proper measure of the computational complexity of such approaches has not been systematically addressed. As systems become larger and more complex, the efficiency of the DEVS simulation algorithms becomes a major issue. Such efficiency can be tackled from at least the following two angles:

• By experimentally comparing implementations using standardized or ad-hoc benchmarking study cases -and preferably large-scale simulation models. [START_REF] Franceschini | A survey of modelling and simulation software frameworks using Discrete Event System Specification[END_REF][START_REF] Van Tendeloo | An evaluation of DEVS simulation tools[END_REF] The results will depend on hardware and software factors such as the processor(s) used, the memory access time, the programming language, or the compiler used. • By theoretically comparing algorithms with each other. [START_REF] Knuth | Selected Papers on Analysis of Algorithms[END_REF][START_REF] Balakirsky | Comparing algorithms: Rules of thumb and an example[END_REF] The advantage of this approach over the preceding one is its independence vis-à-vis the material factors. Its disadvantage is the relative hardness to properly elaborate the theoretical measurement.

This work addresses the theoretical approach to DEVS simulation efficiency. As the DEVS simulation scheme is a tree-based architecture (see Fig. 1), it is important to observe the behavior of the simulation process when the number of simulation components grows width-wise and depth-wise. The term width refers to the number of components per level of the simulation tree, while the term depth refers to the number of levels of such a tree. The paper tries to provide a generic method to such an analysis that can apply to any DEVS abstract simulator. At this stage, it is important to notice that, similar to the complexity analysis of sorting algorithms, the computational complexity measure of DEVS abstract simulators is not intended to provide efficient implementations, but to allow us understand how an algorithm behaves as the input (i.e., width and depth of the simulation tree generated from the structure of the DEVS model to be simulated) grows larger. Therefore, this computational complexity cannot be used for predicting/comparing specific implementations. For the latter purpose, experimental results need to be used (taking also into account hardware/middleware-based optimizations).

In Sec. 2, related works are examined. In Sec. 3, the DEVS simulation architecture is outlined, and the reference model and notation that will be the basis for our proposal are introduced. In Sec. 4, the approach used for determining the computational complexity of DEVS abstract simulators is described. In Sec. 5, the method is applied to various well-known variants of DEVS simulation algorithms. The results obtained are discussed in Sec. 6. Section 7 concludes the paper and gives perspectives for future work.

Related Works

Performance evaluation and measurement constitute a very important aspect of model design. [START_REF] Sun | Performance measurement of dynamic structure DEVS for large scale cellular space models[END_REF] Some research works have addressed the experimental evaluation of DEVS simulation performances; [START_REF] Franceschini | A survey of modelling and simulation software frameworks using Discrete Event System Specification[END_REF][START_REF] Van Tendeloo | An evaluation of DEVS simulation tools[END_REF][START_REF] Risco-Martin | Reconsidering performance of DEVS modelling and simulation environments using the DEVStone benchmark[END_REF][START_REF] Glinsky | DEVSTONE: a benchmarking technique for studying performance of DEVS modelling and simulation environments[END_REF][START_REF] Wainer | Studying performance of DEVS modelling and simulation environments using the DEVStone benchmark[END_REF] few have addressed the theoretical evaluation of their computational complexity.

Among the experimental works, a noticeable one is DEVStone, a benchmark and model generator used for performance evaluation [START_REF] Risco-Martin | Reconsidering performance of DEVS modelling and simulation environments using the DEVStone benchmark[END_REF][START_REF] Gutierrez-Alcaraz | Experiences with the DEVStone benchmark[END_REF] which enables standardized and exhaustive ways to compare different DEVS-based simulation environments. A performance evaluation of different and successive versions of CD++ (a C++ implementation of DEVS-based cellular automata simulation) has been conducted using DEVStone. [START_REF] Balakirsky | Comparing algorithms: Rules of thumb and an example[END_REF][START_REF] Sun | Performance measurement of dynamic structure DEVS for large scale cellular space models[END_REF] In another work, the authors presented a revision of DEVStone [START_REF] Risco-Martin | Reconsidering performance of DEVS modelling and simulation environments using the DEVStone benchmark[END_REF] by introducing new equations and model benchmark, which they used to test five different DEVS-based environments with respect to execution performance and memory footprint on two hardware platforms. In Ref. 13, the authors considered using DEVStone to conduct a technical survey of some of the major DEVS-based simulators and software frameworks, based on a set of criteria to capture the main features of each of them, and furthermore provided a classification upon the tests conducted using these criteria. In Ref. 14, the authors presented a decisional work which enables modelers to choose the best-suited tool to solve their problem and furthermore depicted possible ways to extend (or upgrade) the studied tool.

An early work on theoretical approach to DEVS complexity is given in Ref. 22. For design consideration, a performance evaluation model is considered to find an optimal model decomposition that can be mapped onto a multiprocessor. Another work, [START_REF] Zeigler | Mapping hierarchical discrete event models to multiprocessor systems: Concepts, algorithm, and simulation[END_REF] which improved upon Ref. 22, was considered by exploiting the parallelism natural to discrete models in general; in addition, the work proposed a complexity computational measure that is used to express the speed-up factor gotten from an optimal model decomposition. This work was done on the basis of the classic DEVS (CDEVS) formalism [as opposed to parallel DEVS (PDEVS); see Appendix A]. It is worth mentioning that distributed simulation in Ref. 23 refers to parallel/concurrent execution rather than the proper "distributed execution" as the term "distributed simulation" holds a much broader meaning. That is, a true distributed simulation execution would require the use of synchronization mechanisms, [START_REF] Jefferson | Virtual time[END_REF][START_REF] Fujimoto | Parallel and distribution simulation systems[END_REF][START_REF] Jafer | Conservative vs. optimistic parallel simulation of DEVS and Cell-DEVS: A comparative study[END_REF] namely optimistic and conservative mechanisms, to ensure execution correctness. More recently, simulation-based activity metrics have been used to corroborate the choice of a DEVS simulation model among a set of family models based on analytic-based activity metrics. [START_REF] Capocchi | Discrete-event simulation model generation based on activity metrics[END_REF] The activity notion for a DEVS system is commonly referred to as the number of execution of its transition functions. However, the authors agreed that getting a good activity metric of models (in order to evaluate their performance) after the simulation is very difficult. Therefore, they do not provide any methodological guideline for the general case of any DEVS model.

Some works not theoretically driven cut into a theoretical performance analysis as a preamble of an experimental work. A performance measure, which relates to the number of atomic models, internal transition, and external transition, was proposed in Refs. 18, 19, and 28. It enabled to represent formulas that can be used to obtain a theoretical execution time, which was compared to the experimental execution time to derive the eventual overhead. A comprehensive performance measurement was proposed in Ref. 17 for DSDEVS for large-scale cellular space models. The authors considered modeling layer and simulation layer for performance analysis and proposed quantitative work on performance by considering the overhead inherent to adding or removing model during simulation in the dynamic simulation scheme. To achieve significant speed-up in execution time, the authors in Ref. 29 proposed an extension of the DEVS formalism to exploit the capability of executing hierarchical and modular DEVS models with the properties of high internal and external event parallelism.

DEVS Simulation Architecture

DEVS is a formal Modeling and Simulation (M&S) framework that enables models to be constructed in a hierarchical and modular composition. In the DEVS formalism, there are two basic types of models: atomic model and coupled model (which is a network of atomic models). DEVS considers not only the model behavior but also the model structure in its formalism. Model behavior has to do with model input and output events and states, while the model structure has to do with the model composition (i.e., atomic or coupled components).

A DEVS model can be expressed in either Classic DEVS 1 or Parallel DEVS. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF] The latter represents an improvement of the former by: (1) removing from the coupled model specification the tie-breaking function (Select function), which sequentializes the simulation execution by selecting at each simulation cycle only one component to be active; and (2) introducing in the atomic model specification the confluent function, which takes care of eventual simultaneous events that may occur during simulation.

The DEVS modeling formalism has been abundantly presented in the literature. The interested readers can refer to Appendix A for a short introduction. For further reading on the DEVS, interested readers can refer to a didactic presentation of

DEVS simulation tree

For simulation executions, a hierarchical DEVS model can be mapped onto a simulation tree, as shown in Fig. 1. DEVS formalism explicitly separates the model from its execution engine. The corresponding coordinators and simulators carry out the simulation of, respectively, the coupled and atomic models. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] The simulation process is dictated by the DEVS simulation algorithm used, to simulate the model of interest. However, each of these simulation algorithms has different message concepts, although following the same semantics.

For simulation executions, a hierarchical DEVS model can be mapped onto a simulation tree, as shown in Fig. 1. DEVS formalism explicitly separates the model from its execution engine. The corresponding coordinators and simulators carry out the simulation of, respectively, the coupled and atomic models. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] The simulation process is dictated by the DEVS simulation algorithm used, to simulate the model of interest. However, each of these simulation algorithms has different message concepts, although following the same semantics.

Normalized representation for DEVS simulation algorithms

Various DEVS simulation algorithms exist in the literature. Let us introduce the following definition of our own. Definition 1. A DEVS simulation algorithm is in the normal form if it is depicted in a triggered-by-message form, i.e., as a set of action batches to be performed, with each batch being guarded by the receipt of a given type of message. This normalized representation of the DEVS simulation algorithms mirrors the event-driven programming paradigm, in which the flow of instructions is determined by the triggering of events (i.e., user actions) as presented by Table 1. This form makes it easy to analytically compute the time complexity of the proposed algorithm, as the complexity related to the receipt of each type of message can be computed in isolation and the complexity of the overall algorithm can be derived from the recursion between the receipt of a type of message and the sending of other types of messages or/and the forwarding of the message received. This forms the foundational principles of our approach. Appendix B presents in their normal form the DEVS simulation algorithms to which our method is applied in Sec. 5.

Generic Approach to Complexity Analysis of DEVS Algorithms

Traditionally, the algorithmic complexity of a particular problem is tied to the size of the problem's input (as the length of an array to be sorted). In this work, the input is the DEVS simulation tree, and the size is given by two parameters: width, which corresponds to the number of children per coordinator, and depth, which corresponds to the number of levels of the DEVS simulation tree. Another important aspect is the cost of communication between two computational nodes, in the case the simulation algorithm is distributed over a computer network. The following subsections present the reference architecture and the notation introduced to support our strategy for the determination of the DEVS computational complexity.

Reference architecture

The study of the computational complexity of an arbitrary DEVS model is inherently complex. In order to ease the computation of the time complexity, a balanced DEVS simulation tree is used, where each coordinator (except the root coordinator) possesses the same number N of children, as shown by Fig. 2 (in that case, each coordinator has three children). The use of a balanced tree enables the exhibition of patterns and eases the use of recursion due to the hierarchical nature of the model. Let us assume the simulation tree levels are numbered from 0 to L (= depth), starting from the leaves (therefore, the root coordinator is the L-level node).

Notice that with such an assumption, L is necessarily greater than or equal to 2. For L = 2, the tree is composed of the root coordinator (at level 2), the top-most coordinator (at level 1), and the N simulators (at level 0). Also, N is necessarily greater than or equal to 2, as having a node with only one child does not make sense. To illustrate the principle of recursion, the sequential processing of the part enclosed within the circle at the bottom left of Fig. 2 is detailed by Fig. 3 in the case of the simulation algorithm proposed in Ref. 5. Let us start with the generation of an output event caused by the reception of an @-message by a (j-1)-level node (a simulator in this case). After computing its output (y), the node will send it to its parent coordinator, which in turn processes it and generates new messages to its influences, and/or forwards it to parent (if the latter is part of influences). Each of the influences other than parent will receive a q-message, and will send back a Dmessage when it is done with the processing of that q-message. The circled area of Fig. 3 is a zoom on the computational process triggered by the receipt of a q-message by a j-level node from its (j + 1)-level parent node. The cascade of computations goes on, with each receipt of message triggering a batch of computations, and each batch of computations generating messages to be sent. Therefore, the computational complexity due to the starting @-message will be the sum of the complexity due to the computation of the y output by the receiving node and the one due to the sending of the message to its parent node. This latter one is computed the same way, leading to a recursion in the computation process.

2350013-6

Reference notations

There are two main approaches to the theoretical measurement of algorithm complexity: time (taken by the processor) and space (in memory). Space complexity is much less important nowadays because computers have tremendous memory. This work focuses on time complexity, which relates to a scarcer resource. The idea behind the concept of time complexity of an algorithm is mainly measuring the total time required by an algorithm to run successfully until its completion solving a particular problem. To support our strategy for complexity analysis, a notation is suggested to be adopted with all DEVS simulation algorithms. The following are defined:

M :
Set of all types of messages of the algorithm. τ (m, j) : Time required to completely process a message of type m at levelj of the simulation tree. This may include the time needed to do local computations, send some follow-up messages to children, and receive feedback from them. In some other cases, only a local computation is done (and possibly, feedback is given to parent coordinator). t com : Time required to exchange a message between two components located at contiguous levels (i.e., parent-to-child or child-to-parent message exchange).

Strategy to complexity analysis

A DEVS simulation is a repetition of computation cycles until the simulated duration expires or a stopping condition occurs. The main idea in this work is to consider that the computational cost of each of these cycles is decomposable into the cost of computations resulting in the sending of the message that initiates this cycle by the root coordinator and the cost of computations of the receipt of that message by the top-most coordinator. The latter is in turn decomposable into the cost of computations resulting in the sending of messages to the children of the topmost coordinator and the cost of computations of the receipt of those messages by the top-most coordinator's children. Recursively, this computational approach is applied until a simulation cycle is complete, i.e., when the root coordinator is ready to initiate a new recursive thread of computations.

By considering the time it would take for a message to transit between two adjacent nodes, emphasis is made on the fact that the message-passing process also affects the execution performance, [START_REF] Van Tendeloo | Introduction to parallel DEVS modelling and simulation[END_REF] particularly if the model is distributed over a computing network with high latency.

The normal form provides us with a common reading grid for consistently interpreting all DEVS simulation algorithms, while the reference notation provides a unified notation for a unique computational basis. Given a DEVS algorithm in its normal form, our strategy is the following:

Step 1. For each m M , do calculate τ (m,0) for top-down messages, and τ (m,L) for bottom-up messages, in the worst-case scenario.

Step 2. For each m M , express τ (m,j) in the worst-case scenario either as a function of τ (m,0) for top-down messages or as a function of τ (m,L) for bottom-up messages.

Step 3. Express C(N , L) as a function of N , L, and τ (m,j), m ∈ M .

Step 4. Replace in C(N , L) each τ (m,j) by its corresponding function of τ (m,0) or τ (m,L).

Step 5. The width-wise complexity of the algorithm is given by C(N ,L) as a function of N when N grows infinitely, and the depth-wise complexity of the algorithm is given by C(N ,L) as a function of L when L grows infinitely.

The worst-case scenario is considered here, i.e., at a given simulation cycle all simulators are imminent and influence each other and their parents as well. Other scenarios (such as best-case and average-case scenarios) can also be considered without changing the global scheme of the approach.

Application

Let us illustrate here the application of our approach with the algorithm proposed in Refs. 4 and 5. The five steps previously presented are detailed hereafter.

Step-by-step complexity measure

Step 1. For each m M , do calculate τ (m, 0) for top-down messages, and τ (m, L) for bottom-up messages, in the worst-case scenario:

• τ (@, 0) = 2t com + τ (y, 1) + τ (D, 1) + a, where a includes time to compute an output event;

• τ (y, L) = 0;

• τ (q, 0) = t com + τ (D, 1) + b, where b includes time to manipulate the bag;

• τ ( * , 0) = max(t int , t ext , t con )+t com +τ (D, 1)+c, where c includes time to compute t L and t N ;

• τ (D, L) = t com + τ ( * , L -1); • τ (D, L) = 0
Note that D' is introduced such that (D',t) = (D, tN), as the algorithm distinguishes between the receipt of (D, t) and of (D, tN) by the root coordinator. Also, if the recursion had to continue over the entire simulation run, it would give τ (D, L) = t com + τ (@, L -1), since each receipt of D' by the root coordinator initiates a new cycle by the sending of an @-message to the top-most coordinator.

Step 2. For each m M , express τ (m, j) in the worst-case scenario either as a function of τ (m, 0) for top-down messages, or as a function of τ (m, L) for bottomup messages:

• For m= @, τ (@, j) = N τ(@, j -1) + (N + 1)t com + τ (D, j + 1) + N d,
where d includes time to handle the synchronize set. Applying the same recursive relation, it gives τ (@, j -1) = N τ(@, j -2) + (N + 1)t com + τ (D, j) + N d Therefore, replacing τ (@, j -1) by this formula in the first expression leads to τ (@, j) = N (N τ(@, j -2) + (N + 1)t com + τ (D, j) + N d)

+ (N + 1)t com + τ (D, j + 1) + N d Hence τ (@, j) = N 2 τ (@, j -2) + N (N + 1)t com + N τ(D, j) + N 2 d + (N + 1)t com + τ (D, j + 1) + N d Arranged differently, it gives τ (@, j) = N 2 τ (@, j -2) + (N + 1)(N + 1)t com + N τ(D, j) + τ (D, j + 1) + (N 2 + N )d
Let us replace in this last equation,τ (@, j -2) by the recursive relation initially established but applied at level j-2, i.e., N τ(@, j-3)+(N +1)t com +τ (D, j-1)+N d, which gives

τ (@, j) = N 3 τ (@, j -3) + N 2 τ (N + 1)t com + N 2 (D, j -1) + N 3 d + (N + 1)(N + 1)t com + N τ(D, j) + τ (D, j + 1) + (N 2 + N )d Arranged differently, it gives τ (@, j) = N 3 τ (@, j -3) + (N + 1)(N 2 + N + 1)t com + N 2 τ (D, j -1) + N τ(D, j) + τ (D, j + 1) + (N 3 + N 2 + N )d
The recursive relation that will successively be applied in every last equation obtained is

τ (@, j -k) = N τ(@, j -k -1) + (N + 1)t com + τ (D, j -k + 1) + N d, until j -k = 1, which also means k = j -1. That way, it gives τ (@, j) = N jτ(@, 0) + (N + 1)(N j -1 + N j -2 + • • • + N + 1)t com + N j -1τ (D, 2) + N j -2τ (D, 3) + • • • + N 2 τ (D, j -1) + N τ(D, j) + τ (D, j + 1) + (N j + N j -1 + • • • + N 3 + N 2 + N )d
It is known that

N j-1 + N j-2 + • • • + N + 1 = (N j -1) N -1
Therefore, one can write τ (@, j) = N j τ (@, 0) + (N + 1)

(N j -1) N -1 t com + j k=1 N j-k τ (D, k + 1) + N (N j -1) N -1 d • For m = y, τ (y, j) = N τ(q, j -1) + (N + 1)t com + τ (y, j + 1) + N d + e,
where e includes time to handle the synchronize set and other local variables. Therefore,

τ (y, j + 1) = N τ(q, j) + (N + 1)t com + τ (y, j + 2) + N d + e
Replacing τ (y, j -1) by this expression in the equation giving τ (y, j), it gives τ (y, j) = N τ(q, j -1) + (N + 1)t com + N τ(q, j)

+ (N + 1)t com + τ (y, j + 2) + N d + e + N d + e
Arranged differently, it gives τ (y, j) = N τ(q, j -1) + N τ(q, j) + τ (y, j + 2) + 2(N + 1)t com + 2N d + 2e

Replacing τ (y, j + 2) by N τ(q, j + 1) + (N + 1)t com + τ (y, j + 3) + N d + e, it gives τ (y, j) = N τ(q, j -1) + N τ(q, j) + N τ(q, j + 1) + (N + 1)t com + τ (y, j + 3)

+ N d + e + 2(N + 1)t com + 2N d + 2e
Arranged differently, it gives τ (y, j) = N τ(q, j -1) + N τ(q, j) + N τ(q, j + 1) + τ (y, j + 3)

+ 3(N + 1)t com + 3N d + 3e
Repeating the same recursive replacement until reaching the level L, it leads to τ (y, j) = N τ(q, j -1) + N τ(q, j)

+ • • • + N τ(q, L -2) + (y, L) + (L -j)(N + 1)t com + (L -j)N d + (L -j)e
Therefore, one can write

τ (y, j) = τ (y, L) + N L-j-1 k=0 τ (q, j + k -1) + (L -j)((N + 1)t com + N d + e) • For m = q τ (q, j) = b • For m = * τ ( * , j) = N τ(qj -1) + 2N t com + N τ( * , j -1) + N f,
where f includes time to handle the local variables. By replacing τ ( * , j -1) by N τ(qj -2) + 2N t com + N τ( * , j -2) + N f , it gives τ ( * , j) = N τ(q, j -1) + 2N t com + N 2 τ (q, j -2)

+ 2N 2tcom + N 2 τ ( * , j -2) + N 2 f + N f Arranged differently, it gives τ ( * , j) = N 2 τ ( * , j -2) + N τ(q, j -1) + N 2 τ (q, j -2) + 2N (N + 1)t com + N (N + 1)f
By repeating the recursive replacement, it gives

τ ( * , j) = N jτ( * , 0) + N τ(q, j -1) + N 2 τ (q, j -2) + • • • + N j τ (q, 0) + 2N (N j-1 + • • • + N + 1)t com + N (N j-1 + • • • + N + 1)f As N j-1 + N j-2 + • • • + N + 1 = (N j -1) N -1 , one can write τ ( * , j) = N j τ ( * , 0) + 2N (N j -1) N -1 t com + j k=1 N k τ (q, j -k) + (N j -1) N -1 N f. • For m = D, τ (D, j) = τ (D, j) = tcom + τ (D, j + 1) + N g,
where g includes time to manipulate the local variables. By replacing τ (D, j + 1) by t com + τ (D, j + 2) + N g, it gives

τ (D, j) = t com + t com + τ (D, j + 2) + N g + N g = τ (D, j + 2) + 2t com + 2N g
By repeating the recursive replacement until we reach the level L, it gives

τ (D, j) = τ (D, j) = (L -j)t com + N (L -j)g + τ (D, L).
Step 3. Express C(N L) as a function of N , L, and τ (m, j), m M:

C(N L) = tcom + τ (@, L -1) + h,
where h includes the initialization time (before simulation cycles).

Step 4. Replace in C(N, L) each τ (m, j) by its corresponding function of τ (m, 0) or τ (m, L):

From Step 2, we know (replacing j by L -1) that

τ (@, L) = N L-1 τ (@, 0) + (N + 1) (N L-1 -1) N -1 t com + L-1 k=1 N L-1-k τ (D, k + 1) + N (N L-1 -1) N -1 d Therefore, C(N, L) = tcom + h + N L-1 τ (@, 0) + (N + 1) (N L-1 -1) N -1 t com m + L k=2 N L-k τ (D, k) + N (N L-1 -1) N -1 d.
Replacing τ (@, 0) by t com + τ (y, 1) + τ (D, 1) + a, which we know from Step 1, and each τ (D, k) by (Lk)tcom + N (Lk)g + τ (D, L), which would be from Step 2, gives

C(N, L) = t com + h + N L-1 (2t com + τ (y, 1) + τ (D, 1) + a) + (N + 1)t com + L k=2 N L-k ((L -k)t com + N (L -k)g + τ (D, L)) + N (N L-1 -1) N -1 d
We know from Step 1 thatτ (D, L) = τ (y, L) = 0, and from Step 2 thatτ (D, 1) = (L -1)t com + N (L -1)g + τ (D, L), as well as τ (y, 1) = τ (y, L)

+ N L-2 k=0 τ (q, k) + (L -1)((N + 1)t com + N d + e), so we get τC(N, L) = tcom + h + N L-1 2tcom + N L-2 k=0 (q, k) + (L -1)((N + 1)t com + N d + e) + (L -1)t com + N (L -1)g + a + (N + 1) (N L-1 -1) N -1 t com + L k=2 N L-k ((L -k)t com + N (L -k)g) + N (N L-1 -1) N -1 d
Arranged differently, we get

C(N, L) = tcom 1 + 2N L-1 + N L-1 (L -1)(N + 1) + N L-1 (L -1) + (N + 1) (N L-1 -1) N -1 + L k=2 N L-k (L -k) + h + N L L-2 k=0 τ (q, k) + N L-1 (L -1)(N d + e) + N L((L -1)g + a) + N g L k=2 N L-k (L -k) + N (N L-1 -1) N -1 d
From Step 2, we know that τ (q, k) = b, therefore it gives

C(N L) = t com 1 + 2N L-1 + N L-1 (L -1)(N + 2) + (N + 1) (N L-1 -1) N -1 + L k=2 N L-k (L -k) + N L (L -1)(b + g + d) + N La + N L-1 (L -1)e + N g L k=2 N L-k (L -k) + N (N L-1 -1) N -1 d + h
Although this expression can be organized in different ways, there is not much point in doing a lot more calculations, as what matters is to derive the dominant terms when L (respectively, N) grows infinitely.

Step 5. The width-wise complexity of the algorithm is given by C(N, L)as a function of N when N grows infinitely, and the depth-wise complexity is given by C(N, L) as a function of L when L grows infinitely:

Algorithms of C(N, L)

Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] (tcom

+ α)N L (L -1) + (2tcom + β)N L-1 (L -1) + (2tcom + γ)N L-1 + tcomN + P L k=1 N L-k (L -k)(tcom + δ) + 2tcom + ε Chow 3 N L (L -1)(tcom + α) + N (L-1) (L -1)tcom + N (L-1) (2tcom + β) + N L+1 λ N-1 N (L-1) -1 N-1 -N (L+1) γ N-1 (L -2) + N L ( N L-1 -1 N-1 )δ + N ( N L-1 -1 N-1 )(2tcom + ωε) + tcom + ε Zeigler et al. 2 N L (L -1)(tcom + α) + N (L-1) (L -1)tcom + N (L-1) (tcom + β) + N (L+1) ( tcom +λ N-1 )( N (L-1) -1 N-1 )-N (L+1) ( tcom+γ N-1 )(L -2) + N L ( N (L-1) -1 N-1 )δ + N ( N (L-1) -1 N-1 )tcom + tcom + ε Schwatinski and Pawletta 6 L(L-1)N L-1 (tcom + α) + β(L-1)N L-1 + N L-1 λ + (N L -1) N-1 (2tcom + γ)tcom + δ
When N grows infinitely (L being constant), the dominant terms are the higherdegree monomials in N , as C(N, L) is then a polynomial in N . These terms are

t com N L (L -1) + N L (L -1)(b + g + d) + N La. Therefore,C(N, L) ≈ (t com + α)LN L ,
where α aggregates all costs for manipulating auxiliary variables (such as synchronize set, time variables, etc.).

On the other hand, when L grows infinitely (N being constant), the dominant terms are the higher-degree exponential functions in L. Therefore,C(N, L)

≈ (t com +g)N L L 2 .
Each complexity measure computed is presented by decreasing order first on the power of N , and then on L. The whole idea in arranging terms the way they appear in Table 2 is to quickly identify the most significant terms to be retained when L or N grows. That way, the complexity order can easily be derived from these expressions presented in Table 2.

For example, consider C(N, L) for Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] as established in Sec. 5.1. We can develop it in the following way:

C(N, L) = t com + 2t com N L-1 + N L (L -1)t com + N L-1 (L -1)2t com + N L-1 (N + 1) N -1 t com - (N + 1) N -1 t com + L k=2 N L-k (L -k)t com + N L (L -1)(b + g + d) + N L a + N L-1 (L -1)e + g L-1 k=1 N L-k (L -k) + N L-1 N N -1 d - N N -1 d + h
We can arrange this expression by starting with the terms containing N L (L -1), then terms containing N L-1 (L -1), subsequently the terms containing N L-1 ,and finally the rest. This gives

C(N, L) = N L (L-1)t com + N L (L -1)(b + g + d) + N L-1 (L -1)2t com + N L-1 (L -1)e + 2t com N L-1 + N L-1 N N -1 d + N L-1 (N + 1) N -1 t com - (N + 1) N -1 t com + L k=2 N L-k (L -k)t com + N L a + g L-1 k=1 N L-k (L -k) - N N -1 d + t com + h
Then, some terms can be grouped together to form constant/negligible values (Latin letters).

To derive the complexity order from the resulting expressions (for each algorithm or abstract simulator), only higher-degree monomials in N are to be retained when N grows infinitely (L being constant) since C(N ,L) is then a polynomial in N , and only higher-degree exponential functions in L are to be retained when L grows infinitely (N being constant). Table 3 depicts the complexity for each of the simulation algorithms (i.e., abstract simulators) studied in the work (see the computational details in Appendix C).

Algorithms of Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] and Schwatinski and Pawletta 6 have a better widthwise complexity than those of Chow 3 and Zeigler et al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] For example, for L = 2 (as already mentioned in Sec. 4.1, L is necessarily greater than or equal to 2), the algorithms of Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] as well as Schwatinski and Pawletta 6 have a widthwise complexity of O(N 2 ), while those of Chow 3 as well as Zeigleret al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] have a width-wise complexity of O(N 4 ).

This does not mean that the specific implementations of Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] or Schwatinski and Pawletta 6 necessarily outperform specific implementations of Chow 3 or Zeigler et al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] Rather, it means that if the models to be simulated have more and more immediate sub-components (for example, from a coupled model having 1000 atomic model components to a coupled model having 1000 atomic model components), the degradation of execution performances observed for any specific implementation (whether optimized or not) of Chow 3 or Zeigler et al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] will be worse than the one observed for an implementation (whether optimized or not) of Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] or Schwatinski and Pawletta. [START_REF] Schwatinski | An advanced simulation approach for parallel DEVS with ports[END_REF] Table 3. Complexity order of some DEVS simulation algorithms (L ≥ 2).

Complexity

Chow et al. Width-wise

O(N L ) O(N 2L ) O(N 2L ) O(N L ) Depth-wise O(L 2 N L ) O(LN L ) O(LN L ) O(L 2 N L )
On the contrary, Chow 3 and Zeigler et al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] have a better depth-wise complexity than Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] and Schwatinski and Pawletta. [START_REF] Schwatinski | An advanced simulation approach for parallel DEVS with ports[END_REF] For example, for N = 2 (as already mentioned in Sec. 4.1, N is also necessarily greater than or equal to 2), Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] as well as Schwatinski and Pawletta 6 have a depth-wise complexity of O(L 2 2 L ), while Chow 3 as well as Zeigler et al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF] have a depth-wise complexity of O(L2 L ).

Therefore, if the models to be simulated are more and more decomposed into sub-components and sub-sub-components (for example, from a coupled model having two model components, each of which has two model components in turn, and so on to level 10, to a coupled model having two model components, each of which has two model components in turn, and so on to level 1000), the degradation of execution performances observed for any specific implementation (whether optimized or not) of Chow et al. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF][START_REF] Chow | Abstract simulator for the Parallel DEVS Formalism[END_REF] or Schwatinski and Pawletta 6 will be worse than the one observed for an implementation (whether optimized or not) of Chow 3 or Zeigler et al. [START_REF] Zeigler | Theory of Modeling and Simulation[END_REF]

Experimental Results

Here, our theoretical results are partly assessed with experimental results obtained from a specific implementation of the DEVS abstract simulators. The focus is on Chow et al.'s 4,5 algorithms, which have been implemented in the Java-programmed SimStudio package. [START_REF] Traoré | SimStudio: A next generation modeling and simulation framework[END_REF] Let us consider a disease-spreading situation in a population composed of individuals, each represented by an atomic model. Each individual has a health status, which is either Susceptible, Infected, or Recovered. While a susceptible individual will remain in his/her health status until being infected, an infected individual will recover after five days (incubation time), and a recovered individual will become susceptible after two days (immune time). A susceptible individual gets infected when in touch with two infected individuals.

The DEVS model of such an individual is defined as follows:

M Individual = X, Y, S, δ int , δ ext , δ conf , λ, ta ,
where

• X = {(CON T ACT, v) | v ∈ {Susceptible, Infected, Recovered}},
CONTACT is the name of the unique input port of the model and v is the value received on it;

• Y = {(HEALT H, v) | v ∈ {Susceptible, Infected,
Recovered}}, HEALTH is the name of the unique output port of the model and v is the value sent on it; • S = {(status, σ) | status ∈ {Susceptible, Infected, Recovered}, σ ∈ }, status is the current health status of the individual, while σ is the remaining time to be in that status; • ta(status, σ) = σ;

• δ int (Infected, σ) = (Recovered, 2), δ int (Recovered, σ) = (Susceptible, +∞), time units are in days, and there is no possible internal transition from the Susceptible health status (only a received infection can trigger a transition from that status); • δ ext ((status, σ), e, x) = (Infected, 5) if x contains at least two Infected statuses, (status, σe) otherwise.

Both the width-wise and the depth-wise growths of the model are experimented:

• For the width-wise growth, let us consider L = 2 [for which it is known the widthwise complexity is O(N 2 )], and let us gradually increase the number of individuals in the population. That way, the corresponding simulation tree (which has only a root coordinator at level 2, a top-most coordinator at level 1, and several children at level 0) will gradually have more and more leaves. • For the depth-wise growth, let us consider N = 2 [for which it is known the depth-wise complexity is O(L 2 2 L )], and let us initially decompose the population into two clusters, each of which is made of two individuals. Then after, let us gradually decompose each cluster into two sub-clusters, each of which is made of two individuals. That way, the corresponding simulation tree will gradually have more and more levels, with each node at a given level having two children.

The disease-spreading model is always fully connected, i.e., all components within any given component send their statuses to each other. Figure 4 shows the principle of the evolution of the model's structure, width-wise and depth-wise. Simulation experiments are executed for various scenarios of population's structure, and each experiment lasts for 100 simulated days. Figure 5 depicts the increase of execution time (ET) in milliseconds, as the model grows width-wise (from N =2 to N = 1024), while Fig. 6 compares the experimental results (dotted line) with the theoretical results (straight line) established in the previous sections. The execution performances are averaged over 10 simulation executions for each scenario. Obviously, Fig. 6 assesses the worst-case nature of the complexity measure (as compared to the experimental results), and Fig. 5 assesses whether the experimental performances degrade akin to the theoretical ones.

Figure 7 depicts the increase of execution time as the model grows depth-wise (from L =2 to L = 10) with a fixed value of N (i.e., 2).

Figure 8 compares the experimental and corresponding theoretical results. Note that the model considered in this work is a balanced hierarchical DEVS model, where at each level the number of subordinates per coordinator is the same. In practice, DEVS hierarchical models are not balanced. In such a case, our approach can be used to find the upper and lower bounds to the computational complexity to be determined. That is, for any regular DEVS model, one can find two DEVS balanced models that prove the boundedness of the complexity measure of the regular model in terms of complexity analysis.

Conclusion

This paper presents a generic approach to the determination of the computational complexity of DEVS simulation algorithms, which exploits their recursive nature when they are expressed in a triggered-by-message form. The approach is based on a balanced simulation tree and uses a reference notation that can apply to any DEVS algorithm. The computational load generated by each type of message of the simulation protocol is first computed, and then used to derive the overall load of the algorithm per simulation cycle. The computational complexity results from the asymptotic growth of the width or the depth of the simulation tree. The approach is applied to four DEVS simulation algorithms, and a comparative study of them is provided.

This approach can be adopted with any DEVS simulation algorithm, whether a variant of CDEVS or PDEVS, or even a DEVS extension that is expressed with CDEVS or PDEVS. However, this paper focuses only on sequential simulation, and only on balanced simulation trees. For simulation that employs a nonbalanced DEVS model, the approach provides lower and upper bound estimates of the complexity. For distributed simulation, additional message types have to be introduced in the DEVS algorithm to ensure the necessary synchronization between simulation components that interact remotely. This has no impact on the applicability of our approach, but the number of processors used must also be considered, as well as the way the communication time between two adjacent simulation nodes varies. Examining the case of parallel and distributed DEVS simulation is envisioned as future work. It is worth mentioning that even if the focus of this approach is on finding a way to measuring the computational complexity of a DEVS abstract simulator, rather than improving the existing ones, it can serve as a computational mean to comparing various solutions and demonstrating any improvement from one solution to another. Step 2. For each m ∈ M , express τ (m, j) in the worst-case scenario either as a function of τ (m, 0) for top-down messages, or as a function of τ (m, L) for bottomup messages:

• τ ( * , j) = N τ( * , j -1) + N t com + N c, where c includes time to manipulate local variables,

⇒ τ ( * , j) = N j τ ( * , 0) + N (N j -1) N -1 (tcom + c); • τ (y, j) = N τ(xj -1) + (N + 1)t com + τ (y, j + 1) + N d, where d includes time to manipulate local variables, ⇒ τ (y, j) = τ (y, L) + N L-j-1 k=0 τ (x, j + k -1) + (L -j)((N + 1)t com + N d); • τ (x, j) = N τ(x, j -1) + N t com + N e
, where e includes time to manipulate local variables,

⇒ τ (x, j) = N j τ (x, 0) + N (N j -1) N -1 (t com + e);
• τ (D, j) = f , where f includes time to manipulate local variables.

Step 3. Express C(N, L) as a function of N , L, and τ (m, j), m ∈ M :

τC(N L) == tcom + ( * , L -1) + g,
where g includes time to initialize the simulation (before entering the loop of simulation cycles).

Step 4. Replace in C(N, L) each τ (m, j) by its corresponding function of τ (m, 0) or τ (m, L):

C(N, L) = t com + g + N L-1 τ ( * , 0) + N (N L-1 -1) N -1 (t com + c) ⇒ C(N, L) = t com + g + N L-1 (max(tint,text,tcon) + 2t com + τ (y, 1) + τ (D, 1) + a) + N (N L-1 -1) N -1 (t com +c)
with τ (y, 1) = τ (y, L)+N L-2 k=0 τ (x, mk)+(L -1)((N +1)+N d), and τ (D, 1) = f . Step 5. The width-wise complexity of the algorithm is given by C(N, L)as a function of N when N grows infinitely, and the depth-wise complexity is given by C(N, L) as a function of L when L grows infinitely:

When N grows infinitely (L being constant, and given L ≥ 2),C(N, L) ≈ (λ + δ)N 2L-2 ,andwhenL grows infinitely (N being constant), C(N, L) ≈ (2tcom + α + γ)N L L. Notice the condition L ≥ 2 always holds, otherwise the simulation tree is reduced to the root coordinator and its child (which therefore is necessarily a simulator) and there is no way for N to grow (as there is no possible top-most coordinator).

C.2. Simulation algorithms of Zeigler et al. 2

Step 1. For each m M,do calculate τ(m, 0) for top-down messages, and τ(m, L) for bottom-up messages, in the worst-case scenario:

• τ ( * , 0) = t com + τ (y, 1) + a, where a includes time to compute an output event;

• τ (y, L) = 0; • τ (x, 0) = max(t int , t ext , t con ) + b, where b includes time to manipulate local variables; • τ (D, L) = 0, as each D-message received by the root coordinator initiates a new cycle.

Step 2. For each m ∈ M , express τ (m, j) in the worst-case scenario either as a function of τ (m, 0) for top-down messages, or as a function of τ (m, L) for bottomup messages: Step 1. For each m M , do calculate τ (m, 0) for top-down messages, and τ (m, L) for bottom-up messages, in the worst-case scenario:

• τ ( * , 0) = t int ; • τ (Y, 0) = t com + τ (y, 1) + a, where a includes time to compute an output event;

• τ (x, 0) = t com + τ (D, 1) + max(t ext , t con ) + b, where b includes time to update time variables; • τ (y, L) = 0; • τ (D, L) = 0, as each D-message received by the root coordinator initiates a new cycle.

Step 2. For each m M , express τ (m,j) in the worst-case scenario either as a function of τ (m,0) for top-down messages, or as a function of τ (m,L) for bottom-up messages:

• Step 5. The width-wise complexity of the algorithm is given by C(N, L) as a function of N when N grows infinitely, and the depth-wise complexity is given by C(N ,L) as a function of L when L grows infinitely: When N grows infinitely (L constant),C(N, L) ≈ (L(L -1)t com + L 2 α + L(βα) + 2t com + γβ + λ)N L-1 , and when L grows infinitely (N constant), C(N, L) ≈ (t com + α)L 2 N L .

Fig. 1 .

 1 Fig. 1. (a) Hierarchy of model components and (b) the associated simulators/coordinators.

Fig. 2 .

 2 Fig. 2. Example of a balanced simulation tree.

Fig. 3 .

 3 Fig. 3. Example of sequential processing based on the simulation algorithm defined in Ref. 5.

  t X : Time required for executing a transition function X by the simulator, X = {int, ext, con}. N : Number of children of a coordinator. L : Number of levels of the simulation tree. C(N, L) : Cost of computations generated by the root coordinator in a cycle of simulation (i.e., one of the repetitive threads of computations performed by the root coordinator).

Fig. 4 .

 4 Fig. 4. Experimental model (width-wise and depth-wise growths).

Fig. 5 .

 5 Fig. 5. Execution performances as the model grows width-wise (from N = 2 to N = 1024).

Fig. 6 .

 6 Fig. 6. Theoretical versus experimental results (width-wise growth; from N = 2 to N = 1024).

Fig. 7 .

 7 Fig. 7. Execution performances as the model grows width-wise (from L = 2 to L = 10).

Fig. 8 .

 8 Fig. 8. Theoretical versus experimental result (depth-wise growth; from L = 2 to L = 10).

•

  τ (y, L) = 0; • τ (x, 0) = b where b includes time to manipulate the bag; • τ (D, L) = 0, as each D-message received by the root coordinator initiates a new cycle.

•L-j- 1 k=0τStep 3 .

 13 τ ( * , j) = N t com + N τ( * , j -1) + c, where c includes time to manipulate local variables,⇒ τ ( * , j) = N j τ ( * , 0) + (N j -1) N -1 (N t com + c); • τ (y, j) = N τ(xj -1) + (N + 1)t com + τ (y, j + 1) + N d, where d includes time for variables update, ⇒ τ (y, j) = τ (y, L) + N (x, j + k -1) + (Lj)((N + 1)t com + N d); • τ (x, j) = N τ(x, j -1) + N t com + N e,where eincludes time to manipulate local variables,⇒ τ (x, j) = N j τ (x, 0) + N (N j -1) N -1 (t com +e); • τ (D, j) = N f ,where f includes time to manipulate local variables. Express C(N , L) as a function of N , L, and τ (m, j), m M: τC(N, L) = t com + ( * , L -1) + g, where g includes time to initialize the simulation.

Step 4 . 5 .C. 3 .

 453 Replace in C(N, L) each τ (m, j) by its corresponding function of τ (m, 0) or τ (m, L):C(N, L) = t com + g + N L-1 τ ( * , 0)+ (N L-1 -1) N -1 (N t com + c) ⇒ C(N, L) = t com + g + N L-1 + (max(t int , t ext , t con ) + 2t com + τ (y, 1) + τ (D, 1) + a) + N (N L-1 -1) N -1 (N t com +c), with τ (y, 1) = τ (y, L) + N L-2 k=0 (x, k) + (L -1)((N + 1)t com + N d).Step The width-wise complexity of the algorithm is given by C(N , L) as a function of N when N grows infinitely, and the depth-wise complexity is given by C(N , L) as a function of L when L grows infinitely:When N grows infinitely (L constant, and L ≥ 2),C(N, L) ≈ (tcom + λ + δ)N 2L-2 , and when L grows infinitely (N constant), C(N, L) ≈ (3tcom+α+γ)N L L Simulation algorithms of Schwatinski and Pawletta[START_REF] Schwatinski | An advanced simulation approach for parallel DEVS with ports[END_REF]

Step 3 . 4 . 1 (

 341 τ ( * , j) = N τ(Y, j -1) + N τ(x, j -1) + 2N tcom + c, where c includes time to manipulate local variables;• τ (Y, j) = N τ(Y, j -1) ⇒ τ (Y, j) = N j τ (Y, 0); • τ (x, j) = τ ( * , j) ⇒ τ (x, j) = τ ( * , j) = N j τ ( * , 0) + j k=1 N k τ (Y, jk) + (N j+1 -1) N -1 (2t com + c);• τ (y, j) = t com +τ (y, j +1)+d, where d includes time to manipulate local variables,⇒ τ (y, j) = τ (y, L) + (Lj)(t com + d); • τ (D, j) = t com + τ (D, j + 1) + N e,where e includes time to manipulate local variables, ⇒ τ (D, j) = τ (D, L) + (Lj)(t com + N e). Express C(N , L) as a function of N , L, and τ (m, j), m M: τC(N, L) = tcom + ( * , L -1) + f,where f includes time to initialize the simulation (before entering the loop of simulation cycles).Step Replace in C(N, L) each τ (m, j) by its corresponding function of τ (m, 0) or τ (m, L):C(N, L) = t com + f + N L-1 τ ( * , 0) + L-1 k=1 N k τ (Y, L -1k) + (N L -1) N -1 (2t com + c) ⇒ C(N, L) = t com + f + N L-1 τ ( * , 0) + 2t com + c).

Table 1 .

 1 Pattern of DEVS simulation algorithm description in the normal form.

	Message	Coordinator	Simulator
	Type of message	Computation done by a	Computation done by a simulator
		coordinator at the receipt of a	at the receipt of a message of
		message of this type, and	this type, and sending of new
		sending of new message(s)	message(s)
		or/and forwarding of the
		message received

Table 2 .

 2 Computational loads per cycle of some DEVS simulation algorithms.

Acknowledgments

This work has been partially supported by the African Development Bank (AfDB).

Appendix A. DEVS Modeling Formalism

DEVS formalism was initially known as CDEVS while presenting some limitations to perform parallel implementation. Some of its limitations include a tie-breaking function that handles simultaneously occurring internal transitions of the components of a coupled model, and the fact that it ignores an internal transition function while occurring at the same time with an external input event (collision), in which case the external transition function always takes place. PDEVS has been introduced to alleviate this drawback. [START_REF] Chow | Parallel DEVS: A parallel, hierarchical, modular modelling formalism[END_REF] In this paper, PDEVS is referred to as DEVS.

An atomic DEVS model is defined by the n-tuple X, Y, S, δ int ,δ ext ,δ conf ,λ,ta , where • X, Y , and S are, respectively, the input set, output set, and state set (at any time, the system modeled is in one of the possible states); • ta: S → +∞ 0 is the time advance function (i.e., it gives the lifespan of each state), with +∞ 0 designating the set of nonnegative real numbers, including +∞; • δ int : S → S is the internal transition function (i.e., it is triggered only when the elapsed time in the system's current state s curr has reached ta(s curr ) without the system being disturbed by any receipt of input); • λ :S → Y is the output function (i.e., it computes the output of the system, each time an internal transition is occurring); • δ ext : Q × X → S is the external transition function [i.e., it is triggered only when the system receives an input, while the elapsed time in the system's current state s curr has not reached ta(s curr )], and Q = {(s, e) | s ∈ S, 0 ≤ e < ta(s)} is called the total state; • δ conf : S × X → S is the confluent transition function [i.e., it is triggered only when the system receives an input at exactly the time when the elapsed time in the system's current state s curr has reached ta(s curr )].

The operational semantics of an atomic model is informally described as follows: At the start, the system is in an initial state and remains there until the time specified by ta is exhausted or until the input event is received. In the former case, an internal transition function occurs then the system switches to another state after sending output event as defined by the output function λ. In the latter case, if input event is received before the specified time, then the external transition function is applied. When a collision occurs, i.e., an external event is received concurrently with the elapsed time equal to the time specified by the time advance function, the confluent function is applied in such a way that the system sends output value and changes to a new state.

A coupled DEVS model is a structure:

• X self and Y self are defined the same way as X and Y are for atomic models (self being here a reference to the coupled model, while component models are referred to using the indices such as i, j, or d); For all influencees j of child i do q = Z i,j (y) Send (q, t) to child j Cache j in synchronize set Wait until (D, t) s are received If self is in

Send (y, t) to parent (q, t) Add event q to the bag Add event q to the bag Send (D, t) to parent coordinator (Continued) If tl ≤ t < tN then For all j in I self and all q in bag q = Z self,j (q) Send (q, t) Step 1. For each m ∈ M , do calculate τ (m, 0) for top-down messages, and τ (m, L) for bottom-up messages, in the worst-case scenario:

• τ ( * , 0) = max(t int , t ext , t con ) + 2t com + τ (y, 1) + τ (D, 1) + a, where a includes time to compute an output event and to manipulate local variables;


Tullio Viola

Letícia Carneiro Piton

Alaa Hamadi

Nabiha Chaumeix

Andrea A Comandini

A new generation kinetic model for pyrolytic soot formation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

KEYWORDS: Soot modeling, pyrolysis, shock tube, kinetic mechanism

Despite decades of investigations, soot chemistry still represents a challenge for experimentalists and theoreticians, mainly due to the complexity of the processes involved. The goal of this work is to develop a comprehensive and semi-detailed kinetic model for soot prediction based on an up-to-date PAHs chemistry growth model developed in our serial works on the pyrolysis of a large number of fuels and fuel mixtures using single-pulse shock tube technique at ICARE. The solid-phase sub-model, generated by a Python code, will be based on the sectional method, using pseudo lumped carbon species, called BINs, able to approximate the soot particle distribution by discrete sections of C and H atoms. All the main classes of the soot chemistry will be considered: HACA, inception, growth, coalescence , dehydrogenation and aggregation. The kinetic parameters will be chosen based on reaction classes. The mechanism validation will be done against shock tube laser-based experiments. Corresponding author : tullio.viola@cnrs-orleans.fr

State of the work at present

Detailed profiles

Example of global parameters

  , L. Carneiro Piton, A. Hamadi, N. Chaumeix, A.Comandini CNRS-ICARE, 1C avenue de la Recherche Scientifique 45071 Orléans Cedex 2, France

✓

  Gas phase C1 up to C18 [1] ✓ Soot model based on CRECK model from Politecnico di Milano [2] ✓ Discrete sectional method with pseudo lumped species BINs ✓ Creation of the Soot Mechanism Automated Generator SMAuG framework in Python; ✓ Cantera for tests and comparisons -in progress + PAHs ❑ A new schema for soot class reactions; ❑ Python tools for analysis and managing of kinetic mechanisms; ❑ Validation of the kinetic mechanism on shock tube 28th ICDERS June 19-24, 2022 Napoli, ITALY m = complex refractive index of soot particles 𝑙 = length crossed by the incident ray λ = wavelength of the incident ray I = intensity of the transmitted ray I 0 = intensity of the incident ray Im value = 0.19 NEXT? 𝑑 𝐸𝑆𝐷 = equivalent spherical diameter Simulation @ T=1770K, P=18 bar, C7H8 0.1% in Ar 𝑑 𝐸𝑆𝐷 = 6.37 nm 𝑑 𝐸𝑆𝐷 = 2.02

Acknowledgements

This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 756785). Thanks to Andrea Nobili for help and suggestions on CRECK soot;