Paul Mangold
email: paul.mangold@inria.fr

Alexandre Filiot

Mouhamed Moussa

Vincent Sobanski

Gregoire Ficheur

Paul Andrey

Antoine Lamer

A Decentralized Framework for Biostatistics and Privacy Concerns

Keywords: federated learning, data privacy, biostatistics

Biostatistics and machine learning have been the cornerstone of a variety of recent developments in medicine. In order to gather large enough datasets, it is often necessary to set up multi-centric studies; yet, centralization of measurements can be difficult, either for practical, legal or ethical reasons. As an alternative, federated learning enables leveraging multiple centers' data without actually collating them. While existing works generally require a center to act as a leader and coordinate computations, we propose a fully decentralized framework where each center plays the same role. In this paper, we apply this framework to logistic regression, including confidence intervals computation. We test our algorithm on two distinct clinical datasets split among different centers, and show that it matches results from the centralized framework. In addition, we discuss possible privacy leaks and potential protection mechanisms, paving the way towards further research.

Introduction

The advent of machine learning methods and the ongoing movement towards wide and high-quality data collection have made biostatistics a crucial component in medical research. Constituting large and representative datasets, which are mandatory either to have enough statistical power or to improve models' generalization, is not always feasible within a single medical center. A popular approach is thus to centralize data from multiple centers in one leading site and conduct the study there. With medical data, this centralization is often a practical challenge, as data is sensitive and must be handled within a controlled environment abiding by strong legal and ethical constraints.

An alternative approach, known as federated learning, consists in training statistical models in a decentralized way, leaving the data on each site, running computations locally and communicating aggregated information between centers during the training phase. Such an approach has already been applied to medicine in a few studies, with the goal of preserving the privacy of sensitive data [START_REF] Jochems | Developing and Validating a Survival Prediction Model for NSCLC Patients Through Distributed Learning Across 3 Countries[END_REF][START_REF] Bogowicz | Privacy-preserving distributed learning of radiomics to predict overall survival and HPV status in head and neck cancer[END_REF], as well as data owners' sovereignty.

This work takes a first step towards defining and implementing a decentralized learning framework for medicine, which differs from previous works in that it allows full decentralization, meaning that it does not require any center to play a central role in the computation (although the latter case remains an option). We aim at proving that this framework can produce results virtually identical to the ones obtained in a centralized setting on actual clinical data. To do so, we use two distinct datasets, fit logistic regressions and compute confidence intervals of the estimates. Finally, we put our work into perspective by highlighting some privacy concerns, together with privacy-preserving mechanisms that could address them, depending on desired privacy levels.

Methods

Decentralized Protocol for Logistic Regression

Logistic regression is fit by estimating the parameters that maximize the likelihood over the observed dataset. Iterative algorithms, such as gradient descent, are commonly used to do so. A decentralized version of gradient descent, as described in [START_REF] Lian | Can Decentralized Algorithms Outperform Centralized Algorithms? A Case Study for Decentralized Parallel Stochastic Gradient Descent[END_REF], can thus be used. In this setting, each center runs the following protocol: • Confidence intervals computation. Compute Fisher information on local dataset, and send it to others. Use these values to compute global confidence intervals.

Note that categorical variables are encoded as dummy variables, whose proper encoding requires either a set of agreed-upon values, or extra communications to determine those.

Datasets and Learning Scenarios

The first clinical dataset used in our experiments consists in measurements collected during caesarean sections performed at the Lille University Hospital. We aim at predicting fetal acidosis at birth based on six explanatory variables, including blood pressure drops during the operation. To simulate a multi-centric environment, the 775 records were randomly assigned to four equally-sized (up to one sample) chunks. So as to provide reproducible results, the UCI heart disease dataset [4], available at https://archive.ics.uci.edu/ml/datasets/heart+Disease, was also used. We aim at predicting the presence of a heart disease based on twelve explanatory variables, which mainly encompass clinical measurements at rest and during a controlled physical effort. This data was collected in four distinct medical centers, with variable sample sizes (respectively 303, 261, 130 and 46 records, for a total of 740).

Three different learning scenarios are studied. The "centralized" scenario, in which the entire collated dataset can be used by a single center. The "all alone" scenario, in which each center tries to perform the study using only its local dataset. The "decentralized" scenario, in which centers communicate together without directly exchanging data records, following the protocol detailed in section 2.1.

To compare those three scenarios, we observe the estimated coefficients and their confidence intervals, checking whether they match, and if so, how precisely.

Results

Implementation

Our decentralized framework is implemented as a R package, available under the MIT license at https://gitlab.com/include-project/federate. The developed package handles network communications, and provides a way to simulate decentralized algorithms locally for testing purposes. It currently implements logistic regression with basic gradient descent, but may easily be expanded to comprise new algorithms, as only logical parts need to be re-implemented.

Algorithms are implemented in R with C++ integration using Rcpp. C++ libraries Armadillo and Asio are used for linear algebra and networking, with their respective R bindings RcppArmadillo and AsioHeaders, available at CRAN.

Experimental Results

The three learning scenarios were run on both the caesarean section and heart disease datasets. In the decentralized scenario, the algorithm is run for a few thousands iterations, inducing as many communication rounds. For both datasets, resulting odds ratios and confidence intervals are the same as in the centralized scenario up to 10 -5 precision on each coefficient. Better precision may be achieved at the cost of more communication rounds and tuning. As for the all alone scenario, it results in unsatisfactory estimations in each center, due to insufficient population size. Odds ratios obtained in these various settings are shown in Figure 1. Clinical results on the caesarean dataset match those reported in [START_REF] Kee | Multivariate analysis of factors associated with umbilical arterial pH and standard base excess after Caesarean section under spinal anaesthesia[END_REF]. Scripts for the heart disease dataset are provided in the git repository.

Privacy Concerns

Although decentralized machine learning naturally favors privacy preservation, keeping the data on site does not fully prevent sensitive information leaks. For federated deep learning, [START_REF] Geiping | Near optimal neural network estimator for spectral¨ x-ray photon counting data with pileup[END_REF] show that shared information may reveal parts of the training dataset.

Since logistic regression shares the same underlying optimization procedure as deep learning, it may be vulnerable to similar attacks. This raises major concerns as local datasets are often small (e.g. for studies on rare diseases), the whole purpose being to gather enough data records to achieve statistical significance. Furthermore, individual showing a subset of variables for readability). Circles, triangles and squares respectively represent the "centralized", "all alone" and "decentralized" scenarios. Each color represents odds ratio learned by a center.

Horizontal grey line is 1, and confidence intervals not crossing it suggest a correlation exists between the variable and the outcome.

records are not the only sensitive information that may be revealed: local aggregated values, e.g. mortality rate, may be retrieved, which can expose centers' internal practices.

Discussions

Privacy Improvement Mechanisms

Local aggregates can be protected through secure aggregation [START_REF] Bonawitz | Practical Secure Aggregation for Privacy-Preserving Machine Learning[END_REF], a protocol that consists in adding random masks on sent information that cancel out when computing the result. This yields an exact global average while preventing sent information from being revealed. It could be used during aggregation and confidence interval computations steps. However, this does not protect individual records from leaking. Differentially Private (DP) mechanisms [START_REF] Dwork | The Algorithmic Foundations of Differential Privacy[END_REF] address this problem, by adding noise that blurs individual contribution on shared values, making it almost impossible to guess the presence of an individual in the dataset. This, however, widely impacts results' precision, and obtaining good accuracy while guaranteeing privacy generally requires very fine tuning of algorithms. Such mechanisms could be used at every communication step of our protocol, either before sending values (Local DP) or after their aggregation (Global DP), depending on trusted parties. Table 1 summarizes the impact of these mechanisms and describes who can infer information, thus requiring others' trust.

Conclusion

Our decentralized framework gives very promising results, near-exactly matching those of the centralized scenario when fitting logistic regressions on two distinct clinical datasets. Its design and open-source implementation allow for its re-use, improvement and extension to other learning tasks. We have also identified a set of privacy-preserving mechanisms whose informed use can ease collaborations between clinical data holders.

Figure 1 .

 1 Figure 1. Odds ratios and 95% confidence intervals learned in the three distinct scenarios for the caesarean dataset (left, randomly uniformly split in 4), and heart dataset (right, split across the 4 actual sources,showing a subset of variables for readability). Circles, triangles and squares respectively represent the "centralized", "all alone" and "decentralized" scenarios. Each color represents odds ratio learned by a center.Horizontal grey line is 1, and confidence intervals not crossing it suggest a correlation exists between the variable and the outcome.

Table 1 .

 1 Privacy mechanisms and who can infer what about records and aggregated values from local datasets.Our framework gives accurate results, echoing the conclusions of[START_REF] Bogowicz | Privacy-preserving distributed learning of radiomics to predict overall survival and HPV status in head and neck cancer[END_REF], within an acceptable number of communications rounds. Its modularity allows further experiments with more advanced optimization algorithms. It can also be extended to different learning tasks, including training deep neural networks, e.g. to classify medical images or learn word embeddings from hospital records. Besides, keeping data on site does not guarantee privacy. The latter should thus receive more attention in the future, notably by implementing secure aggregation and differentially private mechanisms. Precisely quantifying required privacy levels is mandatory to make informed choices of protection mechanisms. Therefore, a comprehensive study of effective data leakage appears to be the next step towards this direction. Full decentralization could further improve privacy, by enabling network topologies in which pairs of centers are distanced based on their mutual trust level. It may also lead to developing broader studies, directly leveraging measures from connected devices at patients' homes, or allow learning personalized parameters adjusted to local specificities.

	Mechanism	Who can infer	Data Records	Local Aggregates	Precision
	None	All	Not Protected	Not Protected	Exact
	Sec. Agg.	All	Only Origin Protected	Protected	Exact
	Global DP	Aggregator(s)	Protected	Not Protected	Inexact
	Local DP	No one	Protected	Not Protected	Inexact
	Sec. Agg. + DP	No one	Protected	Protected	Inexact

P. Mangold et al. / A Decentralized Framework for Biostatistics and Privacy Concerns


Nicolas Bus

Fahad Muhammad

Bruno Fies

Ana Roxin

Semantic topological querying for compliance checking

come

INTRODUCTION

In a world where industrial operations are increasingly complex and globalized, there is a growing need for delivering services informing the users about the levels of "safety" and conformity of a given building. While technical experts and audit specialists are moving towards BIM (Building Information Modeling), model-checking software allow to control building operations and processes, starting from the first lifecycle phase -building design. This improves project standard compliance while ensuring related costs remain affordable.

In this study, we will focus on how geometrical and semantic constraints can simultaneously computed in a compliance checking context. This work is part of the French national project "Digitizing building regulations" (« Numérisation des règles ») which aims at converting existing AEC laws and regulations into logical rules.

RELATED WORK

Traditional approaches for model checking

The traditional approach of verification using MVDXML [START_REF]MVDXML specification and support implementation[END_REF] [START_REF] Zhang | Model view checking: automated validation for IFC building models[END_REF] is very limited. Major limitations are identified such as restricted scope of applying conditions and constraints on several branches of an IFC model, poor geometric analysis of an IFC model, lack of mathematical calculations, support of only static verification of a model, etc. Several works introduced a semantic rule engine oriented [START_REF] Fahad | Towards Validation of IFC Models with IfcDoc and SWRL A Comparative Study[END_REF] approach as a viable alternative.

Semantic approach

In the context of BIM (Building Information Modelling), modelling of resources has been identified as an interesting approach for achieving information interoperability (Pauwels, et al., 2011 andFarias, et al., 2014). Existing IFC-related ontologies were conceived as direct syntax mappings between EX-PRESS and OWL languages (Beetz, et al., 2009 andPauwels, et al., 2011). One of the latest and the most solid implementations of an IFC ontology is IfcOWL proposed in (W3C Linked Building Data Community Group, 2014). This version of IfcOWL is also a Candidate Standard (buildingSMART, 2017) for buildingSMART (meaning it is considered as an activity in the process of acquiring international consensus before being submitted to the Standards Committee for a final vote).

In the last years, several initiatives were proposed by the research community, most of them relying on Semantic Web technologies for addressing the issues that appear when relying solely on IFC models. Based on these, buildingSMART International accepted ifcOWL (OWL serialization of the IFC standard) as a candidate standard. Indeed, ifcOWL and its sub-related graph structure allow solving some issues regarding data partitioning, querying and reasoning. By doing so, buildingSMART International (notably its Linked Data Working Group -LDWG) have opened the road to expert modelchecking solutions. Still, at this stage, several questions remain unanswered. Most requirements ad-Semantic topological querying for compliance checking

A. Roxin

Université de Bourgogne, Dijon, France ABSTRACT: Architects and construction engineers needs services for checking their designs against specific standards, regulations and policies. Many works have been done last years to develop checking software. French National project "Digitizing building regulations" aims at formalizing regulation for automatic compliance checking purposes. Our approach is cloud ready, opened, extensible and based on international standard. In this paper we focus on semantic topological aspects to show gains, limits and perpectives of this approach.

dress simultaneously geometrical and semantics aspects.

Taking IFCOWL and IFC to RDF converter as a starting point, we have built transformation services and a methodology that is described in detail in the following chapters.

Collaborative BIM platform

In March 2018, the KROQI platform has been officially launched. It is a collaborative cloud platform for construction SMEs available free of charge to all construction stakeholders including SMEs and VSEs.

It provides document management functionalities, collaborative services, and innovative business services all along the building life cycle.

The MVDXML technology is deployed on the KROQI platform to check the quality of the IFC files that are uploaded by the users. It checks if the structure of the IFC files is compliant with the requirements expressed in the BIM charter.

The semantic checker to be developed in the frame of the «Numérisation des règles » project will complement this first checker. It will be integrated to the KROQI platform and provides extended rule checking capabilities.

METHOD

Validating modeling options

A modeling charter is a document describing how construction project stakeholders agree to model and share information. This document describes data exchange processes, modelling best practices, exports parameters, nomenclatures and classifications to be used. Having all actors of the project share the same modeling rules drastically decrease the number of models' variants. The modeling charter is also used to extend IFC properties with ad-hoc property sets dedicated to specific topics, such as certification level, equipment performances or environmental impact of product. CSTB aims at providing a national modeling charter core to be used as a common starting point when engineering the modeling plan of a project. We are working together with architects, building owner, engineers, audit experts, building quality control agencies and government to reach this goal.

Enriching the model through semantic inferences

The semantic approach allows to describe a building as well as a requirement by using the same atomic fact formalism called triplet. The IFC and so ifcOWL standard vocabulary includes hundreds of classes and properties representing the building physical layer. This vocabulary is not natively de-signed to deal with building performances, requirements and building functionalities.

During the first phase of the "Numérisation des règles" project, a first regulation vocabulary has been extracted from a set of regulatory documents. While trying to align IFC vocabulary with regulatory vocabulary we showed that this alignment cannot be based on a simple bijective relation. Indeed, some regulatory concepts can only be defined with a logical statement aggregating various IFC classes. For instance, in our regulation ontology, a simple regulatory concept like "highest storey" is inferred by comparing the "elevation" (ifcOWL property) values of the all "IfcBuildingStorey" (ifcOWL class) for a given building.

The Regulatory ontology is composed by complex regulatory concepts defined on top of ifcOWL. Regulatory concepts are organized in layers. The ifcOWL vocabulary is the ground layer of a vocabulary pyramid. At the very top of the pyramid (see figure 1), we find the very regulatory-specific vocabulary. A term of a specific layer is defined by using terms belonging to the lower layer. This layered approach is flexible as it keeps high-level definitions simple by increasing the number of terms. The whole process that provides a semantic regulatory view from the IFC source model works as follows:

firstly, we convert IFC model into an ifcOWL triplet database by using the IFC-to-RDF-converter [START_REF] Pauwels | IFC-to-RDF Converter[END_REF] -provided by BuildingSmart-;

A custom class and relations filter is used to remove non-significant information (concepts not involved in any rule). This operation can be compared to a MVD (Model View Definition).

Then a "geometrical preprocessor" renders and computes geometry aspects. This preprocessor is detailed in the following chapter.

Eventually a "semantic preprocessor" infers the model according to the regulatory ontology.

As a major advantage, this last transformation step offers a simplified the graph expressing rela-tions between building objects compared to the native IFC mechanism [START_REF] Pauwels | SimpleBIM: from full ifcOWL graphs to simplified building graphs, Ework and ebusiness in architecture, engineering and construction[END_REF].

For instance, referencing to classification takes at least ten triplets with ifcOWL vocabulary whereas one triplet is sufficient with the regulatory ontology.

Geometrical preprocessing

The regulation describes constraints that involve topological relations between objects such as: is below, is inside, is adjacent to. IFC allows to describe elements with various representation models: extrusion, BREP, boolean operations. To compute topological relations a low-level geometry such as bounding boxes is sufficient.

The building boxes are described by using only 6 semantic relations from the regulation ontology: Xmin, Xmax, Ymin, Ymax, Zmin and Zmax (see figure 2). The geometric preprocessor [START_REF] Preidel | Pre-Processing IFC Building Models for Code Compliance Checking based on Visual Programming[END_REF] [START_REF] Daum | Processing of Topological BIM Queries using Boundary Representation Based Methods[END_REF], infers on geometry to compute triplets materializing bounding box representation of building objects. According to ours tests reasoning on the bounding box representation are reliable and low CPU consuming.

Formalizing regulatory requirements

This section provides details on the methodology used to transform the regulatory constraints, expressed in natural language, into processable queries.

Firstly, regulatory texts are prepared and interpreted by working groups. Experts use a text editor with autocompletion and syntax highlighting capabilities to re-write regulatory texts into semi-formal rules. Each semi-formal rule is supposed to detect non-compliant building elements concerning a very specific aspect. Depending on its complexity, a regulatory rule in text can be divided into several atomic semi-formal rules dealing with complementary aspects. The idea is to keep each semi-formal rule as simple as possible. As a guideline, we suggested that each semi-formal rule begin with "IF" followed by a condition on specific elements and ends with "THEN NON-COMPLIANT".

The second step of this transformation brings the regulatory constraints from the semi-formal to the formal stage. This step is performed by computer scientists familiar with the BIM and with the semantic technologies. The semi-formal stage makes it convenient to build formal rules. Each semi-formal rule is translated to a SPARQL query using the same level of vocabulary coming from the regulatory ontology. The direct use of the regulatory controlled vocabulary within the constraint-queries keeps them understandable for construction experts. Geometrical constraints are described by using the geoSPARQL (OGC standard) relations covering all the possible topological relations between two geometries as introduced by the Trinity College of Dublin [START_REF] Mcglinn | Integrating Ireland's Geospatial Information to Provide Authoritative Building Information Models[END_REF]. Geometry simplification makes it easier to deal with geometry coordinates. This approach could be extended if needed by using functional extensions [START_REF] Zhang | Querying Linked Building Data Using SPARQL with Functional Extensions 7[END_REF].

Checking regulatory requirements as a KROQI service

Automating code-compliance checking consists of chaining conversion algorithm with geometrical and semantic preprocessors introduced in §3.3. Each time new IFC model is submitted to the semantic checker this chain is triggered so that we obtain a triplet database with the right level of details, aligned with the regulatory ontology. The geometrical preprocessor consists of a java code that returns geometrical triplets from IFC. The semantic preprocessor consists of a sequence of semantic forward chaining operations that filters useless information or enhances model with high level vocabulary. Indeed, some IFC terms are almost never addressed by regulatory texts -such as high-level geometry, element history or sensor states so that they can be filtered from the model to reduce the amount of data. Ontology alignment and extensions to IFCOWL are declared in the regulation ontology. Constraintqueries are organized by regulatory topics so that the KROQI user can the set of rules corresponding to his specific needs (fire safety, accessibility, ventilation, acoustics…).

SYSTEM ARCHITECTURE

Technically, the checking service is hosted by a microservice architecture among other BIM service like account manager, document sharing, annotations, instant messaging, BIM viewer. Each service is programmatically independent but can communicate with others through synchronous HTTP/REST queries or asynchronous queue messages.

The checking service is developed according to the SaaS pattern. This service can be executed lo-cally or integrated to any collaborative BIM platform. The service (see Figure 3) can be executed through its own user interface that provides a list of rule set to check.

The REST API published by the service is simple. It consists of a unique method "check" and two parameters. The first parameter defines the path to the IFC or IFCZIP model. Constraint-queries are packed as SPARQL files in a ZIP archive. The archive also contains the metadata for end-users' information purpose. The second parameter defines the path to the ZIP file corresponding to the Constraint-queries set. The result can be returned as XML, BCF, PDF or JSON format according to the accept request-header defined by the REST specification. PDF is used for user display, JSON allow to chain services and BCF can interoperate with BIM-native tools. Whatever the format asked, the result contains for each constraint-query a list of non-compliant building elements identified by their GUID.

The BCF format [START_REF]BCF specification[END_REF] is extremely simple and easy to implement [START_REF] Mahdavi | eWork and eBusiness in Architecture[END_REF]. The basic content is an issue with comments and references to affected objects (using IFC mechanisms for Global Unique ID's). The format also supports status-information, since issues may be discussed by different users during the workflow. The BCF file makes it convenient to display within a standard IFC viewer, on top of a graphical representation of a building element the results of the evaluation. The BCF format is inde-pendent of the IFC schema version, so it can be used with IFC2x3 and IFC4.

RESULTS

We applied this approach on several constraints of the French regulation on both fire safety and accessibility domains. A first version of a French regulation ontology based on previous CSTB research works has been adapted [7]. A group of domain experts had analyzed, interpreted and finally converted a dozen of regulation texts into about one hundred semi-formal constraints implementing concepts of the regulatory ontology vocabulary. The regulation ontology was then extended with experts' new concepts suggestions.The following paragraphs illustrate in detail two specific use cases.

Checking compliance of fire safety regulations

The fire safety domain regulation put some constraint on structure element performances. Indeed, from the fire safety perspective, the structure elements of the building must bear load during at least one hour in case of fire. This constraint differs according to the building height. On a fire safety point of view, the building height is measured from the ground until the floor of the latest storey. Translated as a SPARQL query, this constraint implies various high-level concepts such as: highest level (storey with highest elevation), storey floor (lower slab of a storey), structure element (building element that bear load), loadbearing duration in case of fire. The following sequence is triggered by the constraintquery:

1. The building height according to the fire safety standard is calculated on the fly 2. The material performance threshold is determined according to the building height 3. All structure elements are retrieved and their performances are checked 4. Finally, the query returns all structure elements with fire resistance below the threshold.

Fire safety regulation needs a circulation graph to be checked (see figure 4). This circulation graph can be inferred from the IFC model by using the IfcRelSpaceBoundary relation between spaces and openings. Something some expected relations are missing in the IFC model. To complete information we also infer spaces adjacencies with a fuzzy intersection rule. The fuzzy intersection is a clash detection rule with tolerance on boundary boxes.

This architecture has been implemented on the "French national BIM platform" providing a full automated checking process on various fire safety and accessibility topics.

Checking accessibility requirements

According to French building regulation, a Water Closet (WC) seat must have a free space of 0.8 x 1 meters on its left OR on its right side. The WC is a high level regulatory concept corresponding to an IfcFlowTerminal [START_REF]IFC 2x3, Specification[END_REF] with an IfcFlowTerminalType [START_REF]IFC 2x3, Specification[END_REF] with a predefined type equal to WCSEAT [START_REF]IFC 2x3, Specification[END_REF]. Each free space is classified by using the high level regulatory concept FreeSpace (see figure 5). It represents a virtual object (not physical) with three dimensions (bounding box). According to the regulation, this object must not intersect any other physical element of the building. How these constraints are checked is detailed by the following sequence:

1. For each WC (water closet seat), two FreeSpace elements -and their bounding box-are created on-the-fly. One on the left and one on the right.

2. For each FreeSpace, the query retrieves all building elements that intersect.

3. Finally, the query returns all WC with at least one FreeSpace intersecting at least one physical building element.

In this use case, a bounding box representation of a building element is accurate enough to detect clashes with the WC. Yet, the bounding box must be oriented to be able to locate left and right sides. This orientation can either be explicit or inferred be comparing dimensions.

In addition to the list of elements that break rules, an explanation query displays the list of elements that intersects with a FreeSpace. During our test we manage to detect a handrail that was settled in a wrong location.

Test set

This approach has been validated on various IFC models from the simpliest to the heaviest by using a widely used triplestore implemented in Java with forward and backward chaining capabilities. While the whole process takes around 7 minutes for the Liberty model (2.1 M triplet) on a quad core i5 CPU at 2,5Ghz, it takes 10 minutes for the HITOS model (13.6 M triplets). IFC conversion to RDF is the most consuming phase and should be optimized.

CONCLUSION

We showed that our semantic topological approach allows addressing significant issues raised by the computation of building regulation requirements. We've built, explained and executed queries validating the approach by using portable SPARQL syntax and a triple store implementing W3C recommendations.

Obviously, this approach does not address all constraints defined in the regulations. For instance, when considering a complex corridor with multiples aisles -and their different related widths -the bounding box approximation will output erroneous results while computing intersection between spaces.

DISCUSSION AND FUTURE WORK

Enriching the regulatory ontology

Web semantic technologies bring a powerful formalism to express terminology and constraints. Yet many BIM platform practitioners do not have full understanding of web semantic technics which limit the number of contributors to a regulation ontology. We believe that it is possible to design tools that support visual query building and enable quick, easy and reusable presentation of terminologies and constraints.

Explaining results

Assigning an object to a rule to point out noncompliance can be insufficient to explain complex cases. The EXPLAIN capabilities of SPARQL could be used to give information on contexts and causes, others elements, properties and relation involved in the computation. Ideally the system should give advices on how to solve compliance issues. Topological requirement could be visualized within the model to point out obstructions, clashing or wrong dimensions. We could materialize those requirement as extra elements and then use the BIM-Snippet brought by BCF 2.0 to return those small partial models in a standard way.

In addition, as each BCF-issue can carry a camera/viewport and even a screenshot of the situation those elements could be generated by the checker itself to illustrate the issue.

Expressing compliance on partial elements

For complex elements we note that compliance should be expressed on a single part of the element. For instance one aisle of a complex corridor could be less large than required. As the corridor is generally described as a single space we need a new formalism to assign the issue to a specific part of the space.

Improving topological treatments

Next, we should improve the geometrical preprocessor to compute and reason on topological quantities such as intersection volume, distance between elements or projection surfaces. A full integration of geoSPARQL using WKT description of geometric could help.

Standardizing regulations ontology

To reinforce domain interoperability and build solutions on international and acknowledged standards, we envision to improve the depicted regulations ontology by considering the latest works done in the W3C Linked Building Data Community Group (CG) such as the Building Topology Ontology (BOT), or the ontology defining the core concepts of a building and OntoBREP -ontology for CAD Data and Geometric Constraints.

Figure 1 .

 1 Figure 1. Regulation vocabulary layers

Figure 2 .

 2 Figure 2. bounding box corresponding to the geometry of a chair

Figure 3 .

 3 Figure 3. Checking Service Architecture

Figure 4 .

 4 Figure 4. Circulation graph inferred for IFC model

Figure 5 .

 5 Figure 5. Materialized free space required around WC seat

  N. Bus, F. Muhammad & B.FiesCentre Scientifique et Technique de Bâtiment, Sophia-Antipolis, France

Table 1 .

 1 Models tested

	______________________________________________
	Building name	IFC size ____________	RDF size _____________
	Mb ______________________________________________ Spaces M triplets
	Chanteloup	17	166	0.8
	AC9 Haus	4	26	1.2
	CSTB Building	15	59	2.7
	EM Building	2	17	0.4
	Model 23	29	17	1.9
	Hitos	62	243	13.6
	Liberty _____________________________________________ 12 21 0 2.1

Figure 6. The "Liberty" model displayed in a IFC viewer


Carlotta Donadello

Thi Nhu

Thao Nguyen

Ulrich Razafison

On the mathematical modeling of vole populations spatial dynamics via transport equations on a graph

Keywords: Age-structured population models, Transport equations, Dispersal, Finite-volume schemes, Numerical simulations AMS Classification: 35Q92, 92B05, 92D25, 65M08

In this paper, we propose two models to describe the spatial dynamic of vole populations together with their finite volumes discretization. The models are based on age-structured transport equations set on a graph. The local evolution of the population occurs at the nodes, while transmission between nodes represents spatial dynamics and is a gradual process in the first model, an instantaneous one in the second. The parameters in the reproduction and mortality rates allows to reproduce the characteristics of different landscapes.

Introduction

Small rodents, as voles, are cornerstone species in most temperate ecosystems as their presence is necessary to the survival of a variety of predators. At the same time, rodents can be vectors of serious disease transmissible to humans, and cause significant loss in agricultural production and storage.

Direct observations performed in several different countries suggest that the populations dynamics of voles is characterized by large density fluctuations related to dispersal, i.e. a collective movement in which a significant part of the juveniles in the colony leave to initiate new settlements. A huge literature is devoted to the mathematical modeling of such populations, with no hope to be exhaustive we defer the interested reader to [START_REF] Perthame | Transport equations in biology[END_REF][START_REF] Thieme | Mathematics in Population Biology[END_REF] and the references therein.

Most of the literature only focus on a two variables dynamics (time and age or time and space) but the observations suggest that the evolution depend on the three of them. The papers [START_REF] Coclite | A PDE model for the spatial dynamics of a voles population structured in age[END_REF] and [START_REF] Berthier | Dispersal, landscape and travelling waves in cyclic vole populations[END_REF], and the monograph [START_REF] Webb | Population models structured by age, size, and spatial position[END_REF] use different PDE techniques to deal with three variables models. In the present paper we follow a partially different idea, inspired from [START_REF] Marilleau | Coupling agent-based with equation-based models to study spatially explicit megapopulation dynamics[END_REF]. In [START_REF] Marilleau | Coupling agent-based with equation-based models to study spatially explicit megapopulation dynamics[END_REF], Marilleaux, Lang and Giraudoux reproduced the first colonization dynamics observed from 1998 to 2010 in the Romanche Valley (Hautes Alpes, France) by an hybrid ODE-Multi Agents model structured as follows : the valley area is decomposed into squared cells of side length 100m, so that at a first approximation the evolution of the population inside each of the cells can be described by an ODE (no spatial dynamics). As soon as the total population in a cell reaches a fixed threshold value, which is a fraction of the cell capacity, dispersal occurs, i.e. some of the juvenile voles "leave" the cell, and aggregate into a migrating vole agent with a specific dynamics also depending on the topography and landscape characteristics of the neighbouring cell agents. In the following we refer to this model as MLG.

The advantage of this approach over the purely PDE one is that simulations are fast and their results are visualized in a format which is accessible to a non specialized audience. At the same time, the description of the population dynamics inside cells is relatively rough, and the set of parameters for which computational results match empirical observations is not organized in a continuum, which makes it difficult to predict the behaviour of the model under perturbations.

In the present paper we exploit the idea of decoupling the local population evolution from its spatial dynamics to construct a simple but effective PDE model structured on a graph whose nodes correspond to the cells in the MLG model. Together with the abstract model, we propose its discretization in the form of a suitable finite volumes numerical scheme. Then we run several simulations showing that our model captures many of the features of voles population dynamics as they have been observed in real life. Finally, in view of the construction of an enhanced hybrid model combining a refined PDE description of the evolution in the cells and the Multi Agent approach to represent dispersal, we also consider a model in which the dynamics inside the nodes is splitted from the dynamics on the graph by using different time scales (instantaneous dispersal).

The paper is organized as follows: in Section 2 we describe the models and motivate our assumptions on the basis of observed population characteristics; In Section 3 we present the corresponding finite volumes numerical schemes and validate them by comparison to explicitely computed exact solutions. Section 4 collects some numerical simulations on simple graphs, which show that the models reproduce the qualitative features of the population dynamics, while Section 5 contains a short summary and the perspective directions of this project.

Description of our models

In this section we propose two models to describe the spatial dynamics of a voles population structured in age. Both models write as transport equations on graph.

Voles' sex ratio at birth is 1-1, [START_REF] Cerqueira | Cyclic changes in the population structure and reproductive pattern of the water vole, arvicola terrestris linnaeus, 1758[END_REF][START_REF] Pascal | Essai de typologie de paramètres démographiques et morphologiques de la fraction colonisatrice d'une population de campagnols terrestres (Arvicola terrestris scherman)[END_REF], so that it is enough to model the female population, which is structured into three age classes: young voles, whose age a varies between 0 and A 1 , juveniles, a ∈ [A 1 , A 2 ], and adults, a > A 2 .

The same observations suggest that A 1 ≈ 18 days and A 2 ≈ 56 days, but for simplicity we assume the normalization A 2 = 3A 1 . Since most voles die out of predation, disease or starvation, at a first approximation we could neglect aging when describing the adult population: in particular it is not necessary to consider a decline of fertility over age.

Nevertheless, to avoid non realistic results when performing simulations on a large time horizon (e.g. 4 years), we define A 3 as the maximal age that a vole individual can attain (different values will be taken into account) and we set our model so that individuals of age larger than A 3 , if any, do not participate any more in the evolution of the system.

As numerical values for age transitions, we fix

A 1 = 1, A 2 = 3 and A 3 = 20.
We partition the natural space into squared cells of side length 100 m, V i , i = 1, . . . , n v . We measure the distances between cells centers and, based on observations, we estimate the average time needed by a vole to cover such distances. Then we define a finite weighted graph Γ = (E , V ), whose set of nodes V = {V 1 , . . . ,V n v } represents the cells, while the set of its edges E consists of the triples (V i ,V j , i j ) where V i = V j belong to V (no loops are allowed) and i j > 0 is the travel time between the two vertices. For simplicity we will always assume that Γ is not directed, therefore i j = j i . In general it is possible to assume that the graph is complete, i.e. that there exists an edge joining any given pair of nodes. From the modeling point of view this is compatible with the fact that some individuals (long range dispersers) can be recaptured at observation sites more than 5 Km apart from each other.

Each node V i is endowed with an initial boundary value problem representing the evolution of the density of voles in V i , ρ i , with respect to the time and age variables t and a,

ρ i : R + × R + → R + (t , a) → ρ i (t , a). (2.1)
The total number of voles in V i at time t is given by the integral

Φ i (t ) = A 3 0 ρ i (t , a) da. (2.2)
The spatial dynamics of the population is represented by transmissions of density from one node to the others. The local spatial dynamics, related to foraging activities is neglected. While most individuals spend almost all of their lives in burrows, where they also breed, some of the juveniles have a significant spatial dynamics during dispersals. The mechanism triggering dispersal is still not completely understood and it is most probably related to the interaction of several factors (see for instance [START_REF] Andreassen | New views on how population-intrinsic and communityextrinsic processes interact during the vole population cycles[END_REF][START_REF] Berthier | Dispersal, landscape and travelling waves in cyclic vole populations[END_REF][START_REF] Cerqueira | Cyclic changes in the population structure and reproductive pattern of the water vole, arvicola terrestris linnaeus, 1758[END_REF][START_REF] Delattre | Vole outbreaks in a landscape context: evidence from a six year study of microtus arvalis[END_REF][START_REF] Giraudoux | Population dynamics of fosserial water vole (arvicola terrestris scherman) : a land use and a landscape perspective[END_REF][START_REF] Ims | Spatial synchronization of vole population dynamics by predators birds[END_REF][START_REF] Johnson | Climatic warming disrupts recurrent alpine insect outbreaks[END_REF][START_REF] Saucy | Juvenile dispersal in the vole arvicola terrestris during rainy nights : a preliminary report[END_REF][START_REF] Sherratt | Use of coupled oscillator models to understand synchrony and travelling waves in populations of the field vole microtus agrestis in northern england[END_REF] and references therein). However, for modeling purposes, it is reasonable to link stress conditions to overcrowding, and say that dispersal begins when the total number of individuals in a colony reaches a threshold value, fixed as a fraction of the carrying capacity of the environment, R.

New colonies most often appear in a radius of a few hundreds meters around the original one, but individuals might disperse over more than 5 Km. Topography, in particular relief, drives the dispersal direction, and during this phase voles' mortality is very high.

We define a stochastic matrix K of size n v × n v . Its coefficients k i j fix which fraction of voles leaving the node V i would head toward V j . Since no loops are allowed, we set k i i = 0 for any i ∈ {1, . . . , n v }. At a modeling level we can imagine that k i j depends, for example, on the height or the landscapes in V i and V j , therefore in general k i j = k j i .

From the modeling point of view it makes sense to consider a time dependent matrix K. An example in this direction is given in Section 4.2.2, but for the simplicity of the presentation, we keep K constant for the moment.

Assume that at time t some of the voles initially at V i disperse. Since voles age and die during dispersal but cannot reproduce, the evolution of voles leaving V i at time t and age a and heading toward V j , is described by the ODE

problem    ∂ s u(s) = -d d i s u(s), u(0) = k i j ρ i ( t , a) (2.3)
where d d i s represents mortality during dispersal. The estimated travel time from V i to V j is the weight i j , therefore the density of voles of age a + i j reaching V j at time t + i j is

k i j ρ i ( t , a) exp(-i j d d i s ) = K i j ρ i ( t , a).
Observe that the matrix K = (K i j ) ∈ R n v ×n v is not stochastic, not symmetric and depends on time if and only if K does.

A model with gradual dispersal (GD model)

We always assume that a dispersal begins from the node V i whenever its total population Φ i reaches the threshold R i . In the first model, the dispersal lasts for a time η > 0 after the moment at which Φ i passes again below R i and the departure rate of juveniles during dispersal is c/η, where c > 0 is a parameter to be chosen. Notice that R i is not the capacity of the node V i but a fraction of it, which means that voles are uncomfortable for Φ i = R i but they are not missing resources.

To introduce the initial boundary value problem at the node V i we consider the balance law

∂ t ρ i + ∂ a ρ i = -d i (t , a)ρ i + P i (t , a), (2.4)
where t and a are respectively the time and age variables, d i is the mortality rate at V i , depending on time and age, and P i (t , a) is a source term which accounts for the possible departure/arrival of voles due to dispersal. The exact form of P i is given below.

Together with (2.4) we consider the initial condition

ρ i (0, a) = ρ 0,i (a), (2.5)
and the boundary condition

ρ i (t , 0) = A i (Φ i (t )) A 3 0 b i (t , a)ρ i (t , a) da, (2.6)
where b i is the reproduction rate of the node V i and the coefficient A i (Φ i ) reproduces Allee's effect, see [START_REF] Murray | Mathematical Biology I, An introduction[END_REF]. This means that if the total population of a colony falls below a certain threshold, reproduction rate also diminishes and might become too small to compensate for mortality. The experimental observations suggest that this situation rarely occurs for voles, so that A i ≡ 1 is a reasonable choice in our simulations. However from the numerical point of view no additional difficulty comes from considering non constant A i , see also the Appendix A for a short additional discussion.

For t > 0 given and i ∈ 1, . . . , n v , we call T i (t ) = τ i 1 , . . . , τ i m i (t ) the (possibly empty) set of times τ i r < t at which a dispersal outbreak took place from V i . Then the term P i (t , a) in equation (2.4) takes the form

P i (t , a) = - c η χ [R i ,+∞) ( sup z∈[(t -η) + ,t ] Φ i (z))χ [A 1 ,A 2 ) (a)ρ i (t , a) + c η j =1,...,n v j =i m j (t ) r =1 δ τ j r (t -j i )χ [A 1 ,A 2 ) (a -j i )K j i ρ j (t -j i , a -j i ), (2.7)
where χ [R i ,+∞) and χ [A 1 ,A 2 ) are the characteristic functions of the intervals [R i , +∞) and [A 1 , A 2 ) respectively, t → δ τ (t )

is the Kronecker delta centered at τ and we use the notation x + = max(x, 0).

A model with instantaneous dispersal (ID model)

The model with gradual dispersal introduced above is very solid as it is based on classical PDE's for which a complete well-posedness theory is available. Moreover, the examples in the next section show that it reproduces the main features of the voles' population dynamics. Nevertheless, in order to bring the model to full size, on a graph counting several hundreds of nodes, one should overcome a number of technical difficulties in parallelization and, especially data storage. For this reason and in view of the possible construction of an hybrid PDE-Multi Agent model similar to the MLG, we propose here a second model which leads to a slightly less precise description of the dynamics, but should be easier to implement on large graphs.

Roughly speaking, we modify the previous model so that whenever the total population at a node V i reaches the threshold value R i , the departure of dispersers takes place instantaneously. As a consequence, all the voles who left V i heading for V j will reach the node V j at the same time. These "instantaneous" departures/arrivals are easier to implement in a finite volume scheme and are not too demanding in term of data storage.

In order to write the model for the node V i , we introduce F i = σ i 0 , σ i 1 , . . . , σ i m i , where σ i 0 = 0 and σ i m i = T , be the (possibly reduced to {0, T }) set of all the times between 0 and the final time T at which, either a dispersal outbreak, or an arrival of voles, took place for the node V i . Then between two times σ i K and σ i K +1 , where K = 0, . . . , m i -1, the model can be written as follows

                                 ∂ t ρ i + ∂ a ρ i = -d(t , a)ρ i , ρ i (t , 0) = A (Φ i (t )) A 3 0 b i (t , a)ρ i (t , a) da, ρ i (σ i K +, a) = χ [0,A 1 )∪[A 2 ,A 3 ] (a)ρ i (σ i K -, a) + 1 -χ [R i ,+∞) (Φ i (σ i K -)) χ [A 1 ,A 2 ) (a)ρ i (σ i K -, a) + j =i j =1,...,n v m j -1 r =1 δ σ j r (σ i K -j i )χ [R j ,+∞) (Φ j (σ i K -j i ))χ [A 1 ,A 2 ) (a -j i )K j i ρ j (σ i K -j i , a -j i ).
(2.8)

Reproduction and mortality rates

Voles' reproduction season lasts typically from April to October, but variations are possible due to weather conditions and, of course, depending on the landscape characteristics. The ideal environment for fossorial voles is grassland where both food and shelter are easily available, but they can also settle in woodlands and farmlands. Mortality also depends on landscape characteristics as it is strongly related to the density and the kind of predators, but is everywhere higher during the winter season as food becomes scarce.

Based on [START_REF] Quéré | Les rongeurs de France[END_REF][START_REF] Krebs | Population fluctuations in rodents[END_REF], the basic reproduction rate of a voles' population is 0.

1 (a) = α •    0.15 exp(K (a -A 1 -ε)) if 0 ≤ a ≤ A 1 + ε, 0.3 -0.15 exp(-K (a -A 1 -ε)) if a > A 1 + ε, (2.9)
where we set K = 10 and ε = 0.2. The parameter α ∈ (0, 1] allows to adapt the average rate to specific situations, as we do in Sections 4.1 and 4.2.2.

Observations reported in [START_REF] Krebs | Population fluctuations in rodents[END_REF] suggest that the average mortality rate during the reproduction season varies between 0.04 and 0.6 per 2 weeks, but the mortality rate for juvenile and newborn individuals is 1 to 3 times larger than for adults. Translated to our setting the value of mortality rate is between 0.05 and 0.77 per time unit. These informations lead us to define a mortality rate function of the following form

d 1 (a) = δ •    0.1 -0.025 exp(K (a -A 2 -ε)) if 0 ≤ a ≤ A 2 + ε, 0.05 + 0.025 exp(-K (a -A 2 -ε)) if a > A 2 + ε, (2.10)
where δ > 1 is a parameter to be tuned.The graphs of b 1 and d 1 for α = δ = 1 appear in Figure 1.

To validate our schemes it is enough to run simulations over a short interval of time, t ∈ [0, 1.3], which corresponds to less than one month in real time. In general however we wish to simulate the population dynamics over one or more years, so that we have to take into account the difference between the reproduction season, lasting from April to October, and the cold season where the reproduction stops and the mortality rate is more important, see for instance [START_REF] Giraudoux | Weather influences m. arvalis reproduction but not population dynamics in a 17-year time series[END_REF].

In our context, one year starting from April 1st corresponds to the period [0, T ], where T = 20, and the reproduction season is [0, T * ], where T * = 12. Given κ > 1 we take

K b (t ) =    1 if 0 ≤ t ≤ T * , 0 if T * < t ≤ T,
and and we compute the reproduction and mortality rates as functions of two variables defined on

K d (t ) =    1 if 0 ≤ t ≤ T * , κ if T * < t ≤ T, ( 2
[0, T ] × [0, A 3 ] as b K (t , a) = K b (t )b 1 (a)
and

d K (t , a) = K d (t )d 1 (a).
(2.12)

The values of α, δ and κ will be specified for each of our examples.

To give a first, basic example of the evolution induced from the parameter functions above we consider the evolution internal to one single node over one year computed by the two models. The initial condition at t = 0 consists of 70 adults equidistributed over their age class. We set α = 1, κ = 1.5 and observe that if we choose δ = 1 several dispersals take place in both models starting at time t 8.1. After T * = 12, the total populations decreases and the only voles remaining at t = T are the adult ones. Their total number is larger than the initial condition ( 90 for the GD model and 80 for the ID model), so we say that the colony persists. On the other hand, if the mortality rate is higher, with δ = 5.5, the population declines, no dispersal occurs and the colony disappears before the end of the simulation.

The models coincide in this latter case, while in the first one we can observe several differences. Essentially, some of the juvenile voles do not leave the node during dispersal in the GD model, therefore the total population remains higher and we see more dispersals for the GD model than for the ID one.

Finite volumes approximations of the models

In this section, we propose finite volumes numerical schemes for the models described above and we validate them by comparison with exact solutions, for which analytical expressions exist and can be computed via the classical method of characteristics, [START_REF] Evans | Partial differential equations[END_REF][START_REF] Perthame | Transport equations in biology[END_REF]. As the models rely on transport equations, the schemes are based on upwind schemes, [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF].

Given a finite time horizon T > 0 and a suitable maximal age A 3 , we consider the computational domain [0, T ] × [0, A 3 ] and let ∆t and ∆a be respectively the constant time and age steps. We set N T = T /∆t , N 1 = A 1 /∆a , N 2 = A 2 /∆a and N A = A 3 /∆a , where s stands for the integer part of s. Then for any 1 ≤ h ≤ N A , we introduce the cells

K h = [a h-1/2 , a h+1/2
), the cells centers a h = (h -1/2)∆a and, for 0 ≤ n ≤ N T , we define the times t n = n∆t . For any given i , j ∈ {1, . . . , n v }, i = j , we approximate the edge (V i ,V j , i j ) by (V i ,V j , L i j ∆t ), where L i j = i j /∆t . Analogously, we approximate the age increase during dispersal by L * i j ∆a where

L * i j = i j /∆a . Given a node V i , 1 ≤ i ≤ n v , we denote by b n i ,h , d n i ,h and ρ n i ,h , for 1 ≤ h ≤ N A and 0 ≤ n ≤ N T , the approximation of the averages of b i (t n , •), d i (t n , •) and ρ i (t n , •) on the cell K h , namely b n i ,h = 1 ∆a K h b i (t n , a) da , d n i ,h = 1 ∆a K h d i (t n , a) da, ρ 0 i ,h = 1 ∆a K h ρ i ,0 (a) da and ρ n i ,h 1 ∆a K h ρ(t n , a) da if n > 0. (3.13)
For each i ∈ {1, . . . , n v } and n ∈ {0, . . . , N T } the total population of voles in V i at time t n is

Φ n i = ∆a N A h=1 ρ n i ,h .
The numerical scheme for the GD model For any given value of the parameter η > 0, minimal time span of any dispersal in our model, we introduce the number of time steps N η = η/∆t and we fix to N η ∆t the minimal duration of a dispersal in the numerical scheme. For n ∈ N fixed and

τ i K in T i (t n ) = τ i 1 , . . . , τ i m i (t n )
, the set of outbreak times taking place from V i before t n , we introduce N τ i K = τ i K /∆t . We next define

Φ n i = sup (n-N η ) + ≤m≤n Φ m i , N i j K = N τ i K + L i j and δ n N i j K =    1 if n = N i j K , 0 otherwise. We finally set d n+1 i ,h = 1 + ∆t d n+1 i ,h + c η χ [R i ,+∞) Φ n i χ [N 1 ,N 1 +N 2 ) (h) ,
which comes from the implicit treatment of the mortality term in (2.4) and the semi-implicit treatment of the dispersal term in (2.7). Then, writing a standard upwind scheme for the transport part of (2.4) and using an explicit treatment of the last term of (2.7) we obtain the following scheme

                             ρ n+1 i ,h = 1 d n+1 i ,h     1 - ∆t ∆a ρ n i ,h + ∆t ∆a ρ n i ,h-1 + c ∆t η j =i j =1,...,n v m j (t n ) K =1 δ n N j i K χ [N 1 ,N 1 +N 2 ] (h -L * j i )K j i ρ n-L j i j ,h-L * j i     , for 1 ≤ i ≤ n v , 1 ≤ h ≤ N A , 0 ≤ n ≤ N T -1, ρ n i ,0 = A i (Φ n i )∆a N A h=1 b n i ,h ρ n i ,h , 1 ≤ i ≤ n v , 0 ≤ n ≤ N T -1.
(3.14)

The numerical scheme for the ID model For n ∈ N, let σ i K be in the set F i = σ i 1 , . . . , σ i m i of times at which either a dispersal or an arrival took place in the node V i before t n . We introduce N σ i K = σ i K /∆t . Writing again a standard upwind scheme for (2.8), we get for each i ∈ {1, . . . , n v },

                                   ρ n+1 i ,h = 1 1 + ∆t d n+1 i ,h 1 - ∆t ∆a ρ n i ,h + ∆t ∆a ρ n i ,h-1 , 1 ≤ h ≤ N A , N σ i K ≤ n < N σ i K +1
,

ρ n i ,0 = A i (Φ n i )∆a N A h=1 b n i ,h ρ n i ,h , N σ i K ≤ n < N σ i K +1 ρ N σ i K i ,h = ρ N σ i K - i ,h -χ [R i ,+∞) (Φ N σ i K - i )χ [N 1 ,N 1 +N 2 ) (h)ρ N σ i K - i ,h + j =i j =1,...,n v m j -1 r =1 δ N σ i K -L j i N σ j r χ [R j ,+∞) (Φ N σ i K -L j i j )χ [N 1 ,N 1 +N 2 ) (h -L * j i )K j i ρ N σ i K -L j i j ,h-L * j i , 1 ≤ h ≤ N A .
(3.15)

Validation

In this part, we validate the implementation of the numerical schemes described above by comparison to exact solutions for which an analytical expression exists.

We consider a simple graph consisting of two nodes V 1 and V 2 and, denoting by L the matrix of the travel times and K the distribution matrix, we set L = 0 0.5 0.5 0 and

K = 0 1 1 0 , (3.16)
which means that, during a dispersal in the node V 1 , all the juvenile voles will move to the node V 2 and vice versa. The initial conditions at the nodes are the following

• At V 1 , we consider 50 juveniles and 170 adults. Individuals are equidistributed over their age class.

ρ 1 (0, a) = 50 A 2 -A 1 χ [A 1 ,A 2 ) (a) + 170 A 3 -A 2 χ [A 2 ,A 3 ] (a). (3.17)
• At V 2 , we consider 120 adults. Individuals are equidistributed over this age class.

ρ 2 (0, a) = 120 A 3 -A 2 χ [A 2 ,A 3 ] (a). (3.18)
For simplicity we assume all the other parameters of our models to be piecewise constant and identical on the two nodes

• Reproduction rate: b(a) = bχ [A 1 ,A 3 ] (a) = 0.3χ [A 1 ,A 3 ] (a) ; • Mortality rate: d(a) = d b j χ [0,A 2 ) (a) + d ad χ [A 2 ,A 3 ] (a) = 0.2χ [0,A 2 ) (a) + 0.1χ [A 2 ,A 3 ] (a) ;
• Dispersal threshold: R = 200 ;

• Mortality rate during dispersal:

d d i s (a) = 0.5.
For the GD model we set c = 1 and we also have to fix the value of the minimal time span of dispersal and we set η = 0.25.

The exact solution for the GD model

In this section we provide a detailed description of the explicit solution in the time interval [0, 1.3] for the GD model with the settings and initial conditions introduced in Section 3.1. The travel time between the two nodes is = 0.5.

Solution at V 1 during the first dispersal, t ∈ [0, t 1 ]
. The given initial condition makes the total population in V 1 at t = 0 larger than the dispersal threshold R and triggers a first dispersal. The density of voles in V 1 is a function of time and age, defined piecewise with respect to age

ρ 1 (t , a) =                          babies 1 (t , a) = be -d b j a A 3 A 1 ρ 1 (t -a, c) dc, a ∈ [0, t ), 0, a ∈ [t , A 1 + t ), juv 1 (t , a) = 50 A 2 -A 1 e -(d b j +1/η)t = 25e -4.2t , a ∈ [A 1 + t , A 2 ), newad 1 (t , a) = 50 A 2 -A 1 e -(d b j +1/η)(t -(a-A 2 ))-d ad (a-A 2 ) = 25e -4.2(t -(a-A 2 ))-0.1(a-A 2 ) , a ∈ [A 2 , A 2 + t ), adults 1 (t , a) = 170 A 3 -A 2 e -d ad t = 10e -0.1t , a ∈ [A 2 + t , A 3 ]. (3.19)
We observe that the total density of juvenile voles of age a who have left V 1 at time t is given by 50

A 2 -A 1 e -d b j t -juv 1 (t , a),
so that the density of voles of age

a leaving V 1 at time t is disp(t , a) = 50 A 2 -A 1 e -d b j t -d b j + d b j + 1 η e -t /η . ( 3.20)
The value of the total population in V 1 , Φ 1 , equals the threshold R at t 0 0.18755. From this moment the dispersal continues for a time η more, and we denote by t 1 = t 0 + η 0.43755 the instant at which it ends.

Solution at V 1 between two dispersals, t ∈ [t 1 , t 2 ].
The function ρ 1 is still piecewise defined with respect to the age variable, but its components are different. The evolution of the component juv 1 is juv 12 . The new components correspond to the babies born at t > t 1 , babies 2 , and to the voles passing from juveniles to adults at t > t 1 , newad 2 .

ρ 1 (t , a) =                  babies 2 (t , a) = be -d b j a A 3 A 1 ρ 1 (t -a, c) dc, a ∈ [0, t -t 1 ), juv 12 (t , a) = juv 1 (t 1 , a -(t -t 1 ))e -d b j (t -t 1 ) , a ∈ [A 1 + t , A 2 ), newad 2 (t , a) = juv 1 (t 1 , a -(t -t 1 ))e -d b j (t -t 1 -(a-A 2 ))-d ad (a-A 2 ) , a ∈ [A 2 , A 2 + (t -t 1 )), ρ 1 (t 1 , a -(t -t 1
))e -d(a)(t -t 1 ) , otherwise.

(3.21)

The total population at V 1 , Φ 1 , increases up to reach again the value R at t = t 2 0.84818. This triggers a second dispersal.

Solution at V 1 for t ∈ [t 2 , A 1 ]
. The second dispersal modifies the dynamics of the juveniles, so that the functions giving the density of newborns and the density of new adults also need to be updated. The evolution of the component juv 12 is juv 123 .

ρ 1 (t , a) =                  babies 3 (t , a) = be -d b j a A 3 A 1 ρ 1 (t -a, c) dc, a ∈ [0, t -t 2 ), juv 123 (t , a) = juv 12 (t 2 , a -(t -t 2 ))e -(d b j +1/η)(t -t 2 ) , a ∈ [A 1 + t , A 2 ), newad 3 (t , a) = juv 12 (t 2 , a -(t -t 2 ))e -(d b j +1/η)(t -t 2 -(a-A 2 ))-d ad (a-A 2 ) , a ∈ [A 2 , A 2 + (t -t 2 )), ρ 1 (t 2 , a -(t -t 2
))e -d(a)(t -t 2 ) , otherwise.

(3.22)

Observe the total population in V 1 reaches again the value R at time t 3 0.94675. This does not produce an immediate change in the dynamics because the difference t 3 -t 2 is less that η, but tells us that the second dispersal will last up to time t 4 = t 3 + η 1.44675.

Solution at V 1 for t ∈ [A 1 , 1.3].
For t > A 1 we have to consider two additional age classes : some of the babies becomes juvenile and start to reproduce, we denote them by newjuv 4 , and babies born after time t = A 1 , babies 4 .

ρ 1 (t , a) =                  babies 4 (t , a) = be -d b j a A 3 A 1 ρ 1 (t -a, c) dc, a ∈ [0, t -A 1 ), babies 1 (t , a), a ∈ [t -t 1 , A 1 ), newjuv 4 (t , a) = babies 1 (A 1 , a -(t -A 1 ))e -d b j (t -a)-(d b j +1/η)(a-A 1 ) a ∈ [A 1 , t ), ρ 1 (A 1 , a -(t -A 1
))e -d(a)(t -A 1 ) , otherwise.

(3.23)

Solution at V 2 before the arrival of voles from V 1 , t ∈ [0, ]. At time t = 0 the total population at V 2 is 120. We will see that no dispersal takes place from this node for t = [0, 1.3]. We use a different font to distinguish the age groups appearing in the solution at node V 2 , while the indices are given in the same way.

ρ 2 (t , a) =            babies 1 (t , a) = be -d b j a A 3 A 1 ρ 2 (t -a, c) dc, a ∈ [0, t ), 0 a ∈ [t , A 2 + t ), adults 1 (t , a) = 120 A 3 -A 2 e -d ad t = 120 17 e -0.1t , a ∈ [A 2 + t , A 3 ]. (3.24)
Solution at V 2 during the arrival of the first group of dispersers, t ∈ [ , + t 1 ]. We have to take into account several new groups of individuals. The new components correspond to the babies born at t > , babies 2 , and to the voles arriving from the node V 1 . Some of them arrive as juveniles, juv 2 , some as adults, adults 2 . Some of the juveniles become adults after their arrival, we call them newad 2 . The voles whose age is in

(A 2 + , A 2 + t ) originated from V 1 ,
but when they arrived they were jounger than A 2 + , we call this group adultsevol 2 . We compute their density using (3.20) and the fact that while travelling from V 1 to V 2 the dispersers solve the ODE y = -d d i s y.

ρ 2 (t , a) =                                      babies 2 (t , a) = be -d b j a A 3 A 1 ρ 2 (t -a, c) dc, a ∈ [0, t -), juv 2 (t , a) = t e -d b j (t -s) e -d d i s disp(s -, a -(t -s + )) ds, a ∈ [A 1 + t , A 2 ), newad 2 (t , a) + adults 2 (t , a) = juv 2 (t -(a -A 2 ), A 2 )e -d ad (a-A 2 ) + t t -(a-A 2 ) e -d ad (t -s) e -d d i s disp(s -, a -(t -s + )) ds, a ∈ [A 2 , A 2 + (t -)), adults 2 (t , a) = t e -d ad (t -s) e -d d i s disp(s -, a -(t -s + )) ds, a ∈ [A 2 + (t -), A 2 + ), adultsevol 2 (t , a) = adults 2 (t -(a -(A 2 + )), A 2 + )e -d ad (a-(A 2 + )) , a ∈ [A 2 + , A 2 + t ), ρ 2 ( , a -(t -
))e -d(a)(t -) , otherwise.

(3.25)

Solution at V 2 for t ∈ [ + t 1 , 1.3].
In the interval of time [ + t 1 , + t 2 ] no more voles arrive from V 1 and the total population of V 2 does not attain the threshold value R. For our validation we only compute the solution up to time

t = 1.3 < + t 2 .
The age groups already present in the dynamics just evolve following the prescribed mortality rates, until time t = A 1 when some on the individuals in babies 1 become juvenile. Therefore for t ∈ [ + t 1 , A 1 ] the function ρ 2 takes the form

ρ 2 (t , a) =            babies 3 (t , a) = be -d b j a A 3 A 1 ρ 2 (t -a, c) dc, a ∈ [0, t -( + t 1 )), newad 3 (t , a) = juv 2 ( + t 1 , a -(t --t 1 ))e -d b j (t -( +t 1 )-(a-A 2 ))-d ad (a-A 2 ) , a ∈ [A 2 , A 2 + t -( + t 1 )), ρ 2 ( + t 1 , a -(t --t 1
))e -d(a)(t -( +t 1 )) , otherwise.

(3.26)

While for t ∈ [A 1 , 1.3] ρ 2 becomes ρ 2 (t , a) =                          babies 4 (t , a) = be -d b j a A 3 A 1 ρ 2 (t -a, c) dc, a ∈ [0, t -A 1 ), babies 1 (t , a), a ∈ [t -A 1 , A 1 ), newjuv 4 (t , a) = babies 1 (A 1 , a -(t -A 1 ))e -d b j (t -A 1 ) , a ∈ [A 1 , t ), newad 4 (t , a) = juv 2 (A 1 , a -(t -A 1 ))e -d b j (t -A 1 -(a-A 2 ))-d ad (a-A 2 ) , a ∈ [A 2 , A 2 + t -A 1 ), ρ 2 (A 1 , a -(t -A 1
))e -d(a)(t -A 1 ) , otherwise.

(3.27)

Figure 3 shows the comparison between the exact and the numerical solutions at the two nodes for the GD model at the final time T = 1.3. We observe a perfect matching between the two solutions where the numerical one was computed with ∆a = ∆t = 0.00625.

The exact solution for the ID model

In this section we provide a detailed description of the explicit solution in the time interval [0, 1.3] for the ID model with the settings and initial conditions introduced in Section 3.1. The travel time between the two nodes is = 0.5.

Solution at V 1 for t ∈ [0, A 1 ]
. The given initial condition makes the total population in V 1 at t = 0 higher than the dispersal threshold R, therefore the juvenile individuals immediately leave the node. The density of voles in V 1 is a function of time and age, defined piecewise with respect to age

ρ 1 (t , a) =            babies 1 (t , a) = be -d b j a A 3 A 1 ρ 1 (t -a, c) dc, a ∈ [0, t ), 0, a ∈ [t , A 2 + t ), adults 1 (t , a) = 170 A 3 -A 2 e -d ad t = 10e -0.1t , a ∈ [A 2 + t , A 3 ].
(3.28)

Solution at V 1 for t ∈ [A 1 , 1.3].
Starting from time t = A 1 a fraction of the babies will become juveniles and their reproductive process has to be taken into account. This leads to the formation of two more age groups: the babies born after time t = A 1 , babies 2 , and the individuals passing from the babies to juveniles, newjuv 2 .

ρ 1 (t , a) =                          babies 2 (t , a) = be -d b j a A 3 A 1 ρ 1 (t -a, c) dc, a ∈ [0, t -A 1 ), babies 1 (t , a), a ∈ [t -A 1 , A 1 ), newjuv 2 (t , a) = babies 1 (A 1 , a -(t -A 1 ))e -d b j (t -A 1 ) , a ∈ [A 1 , t ), 0, a ∈ [t , A 2 + t ), adults 1 (t , a) = adults 1 (A 1 , a -(t -A 1 ))e -d ad (t -A 1 ) , a ∈ [A 2 + t , A 3 ].

Solution at

V 2 for t ∈ [ , A 1 ]. The density of voles reaching V 2 at t = = 0.5 are 50 A 2 -A 1 χ [A 1 + ,A 2 + ] (a)e -d d i s
, then some of them are juveniles, juv 2 , some are adults, adults 2 . Some of the juveniles become adults after their arrival, we call them newad 2 . We call babies 2 , the babies born for t > .

ρ 2 (t , a) =                          babies 2 (t , a) = be -d b j a A 3 A 1 ρ 2 (t -a, c) dc, a ∈ [0, t -), juv 2 (t , a) = 50 A 2 -A 1 e -d d i s -d b j (t -) , a ∈ 1 + t , A 2 ), newad 2 (t , a) = juv 2 (t -(a -A 2 ), A 2 )e -d ad (a-A 2 ) , a ∈ [A 2 , A 2 + (t -)), adults 2 (t , a) = 50 A 2 -A 1 e -d d i s -d ad (t -) , a ∈ [A 2 + (t -), A 2 + t ), ρ 2 ( , a -(t -))e -d(a)(t -) , otherwise.
(3.30)

Solution at V 2 for t ∈ [A 1 , 1.3].
Starting from time t = A 1 a fraction of the babies will become juveniles and their reproductive process has to be taken into account. This leads to the formation of two more age groups : the babies born after time t = A 1 , babies 3 , and the individuals passing from the babies to juveniles, newjuv 3 .

ρ 2 (t , a) =                  babies 3 (t , a) = be -d b j a A 3 A 1 ρ 2 (t -a, c) d c, a ∈ [0, t -A 1 ), babies 1 (t , a), a ∈ [t -, A 1 ), newjuv 3 (t , a) = babies 1 (A 1 , a -(t -A 1 ))e -d b j (t -A 1 ) , a ∈ [A 1 , t ), ρ 2 (A 1 , a -(t -A 1
))e -d(a)(t -A 1 ) , otherwise.

(3.31)

Figure 4 illustrates the comparison between the explicit solution and its numerical solution at the two nodes for the ID model at the final time T = 1.3. We again observe a perfect agreement between the two solutions. The numerical solutions was computed with ∆a = ∆t = 0.00625.

Numerical convergence

We introduce here the following relative L 1 -discrete error for the density in the node

V i , i = 1, 2, at a given time t n e n i = 1≤h≤N a |ρ i (t n , a h ) -ρ n i ,h | 1≤h≤N a |ρ i (t n , a h )| ,
where ρ(t n , a h ) is an exact solution evaluated at the point a h and at time t n .

In tables 1 and 2, we have reported the relative L 1 -discrete errors for both models in the two nodes V 1 and V 2 , computed for different values of age and times steps at two times T = 0.6 and T = 1.2. The results suggest the convergence of the numerical schemes. Moreover, figures 5 and 6 suggest that the order of convergence of the scheme for the GD model is between 0.71 and 1, and the one for the ID model is 1. We can also observe that the loss of accuracy of the scheme for the GD model is mainly due to the approximation of the arrival term in (2.7). Table 2: L 1 -discrete relative errors for the ID model

-6 -5 -4 -3 -2 -6 -5 -4 -3 -2 -1 ∆a e NT 1 e NT 1 Slope 1 (a) V 1 at T = 0.6 -6 -5 -4 -3 -2 -6 -5 -4 -3 -2 -

A first comparison between the two models

In this section we consider once again the initial conditions (3.17 We can observe that the profile of the solution is, at a first approximation, similar to the solution obtained from the same initial condition by the ID model, see Figure 4. In particular we observe that the two dispersals occurring in the fist cell for t > 1 concerne a so small number of juveniles that none of them reaches the second cell, so that the dynamics of the second cell is almost identical in the two simulations.

This means that, even if the ID model is necessarily less precise, it suffices to capture the dynamics we would observe in the GD model for small values of η.

Numerical simulations

In this section, we perform numerical simulations on the two models, in order to show their ability to reproduce the qualitative behavior of vole populations and in particular their spatial dynamics. In this paper we limit our attention to simple toy models as the implementation of the schemes on graph representing large regions would require a large effort in optimization and in parallelization of the codes, which we defer to future investigation.

Everywhere in the following we assume, in the GD model, that c = 1 and that the minimal duration of a dispersal is η = 7/18 that corresponds to one week in real time.

Two-nodes dynamics

In this section, we investigate the evolution of a population of voles on a graph consisting of two nodes assumed to represent two different types of landscapes and the edge linking them. The landscapes we consider are described below:

• The first node represents a crop field, which is a source of abundant food and shelter for voles during summer and spring. This traduces into high reproduction rate and low mortality rate in the warm season. On the contrary mortality becomes high at harvesting time and remains so up to the end of one year cycle.

To reproduce such properties, we consider the reproduction and mortality rate functions defined by (2.9), (2.10) and (2.12) where we assume α = 1, δ = 1 and κ = 2.6.

• The second node represents a forest, where the population of voles is almost stable over several years. The reproduction rate is lower than in the crop field, but the mortality varies less dramatically from warm to cold season.

To reproduce such properties, we consider the reproduction and mortality rate functions defined by (2.9), (2.10) and (2.12) where we assume α = 0.656, δ = 1 and κ = 1.5. We assume that the travel time from one node to the other is = 0. Both models reproduce the cyclic behavior of the population. A characteristic feature of the ID model is the appearance of sharp peaks, as in Figure 8(d). This is due to the fact that voles coming from the node representing the crop field reach their destination at moments where the total population in the node representing forest is close to the dispersal threshold. This sudden population increase, sufficient to trigger dispersal, lasts only one time step. On the contrary, no peaks arise in Figure 8(c) as in the GD model, dispersals and arrivals happen progressively. Moreover, since not all of the juveniles leave, the variation in the total population is lower, see Figures 8(a) and 8(c).

Three-nodes dynamics

In this section we run our models on a complete graph with three nodes, V i , i = 1, 2, 3. This allows us to use a more interesting distribution matrix to reproduce, for example, the fact that voles prefer to move downhill rather than uphill, and that they are more likely to stop at places where resources are more abundant. Since the availability of food may change in time, we consider a time dependent distribution matrix in our second example.

The simulations of this section are performed with ∆a = ∆t = 0.0125.

Scenario 1: Three nodes located at different heights

We assume here that the node V 1 is located at the highest spot and V 3 at the lowest one. For all the nodes, we use the mortality and the reproduction rates defined by (2.9), (2.10) and (2.12) with α = δ = 1 and κ = 1.5. As initial conditions we take φ 1 (0) = 100, φ 2 (0) = 50 equidistributed over the adult age class, and φ 3 (0) = 0, that means that the third node is initially empty. It has been observed that when dispersal occur, voles move more frequently downhill (2 times out of 3), see for instance [START_REF] Delattre | Vole outbreaks in a landscape context: evidence from a six year study of microtus arvalis[END_REF][START_REF] Marilleau | Coupling agent-based with equation-based models to study spatially explicit megapopulation dynamics[END_REF]. In other to take into account this characteristic, we consider the following distribution matrix

K =     0 1 0 1/3 0 2/3 0 1 0     .
In particular, this means that all the dispersals taking place from V 1 or V 3 head toward V 2 , while 2 out of 3 voles leaving V 2 would go to V 3 . We fix the travel times matrix as Figures 9 show the evolution of the population for both models over three years. We can observe numerous outbreaks from V 2 , which are due to the fact that all the voles that disperse from V 1 and V 3 arrive in this node. Moreover from Figure 9(c), we can see that, due to the relatively long duration of dispersal, the total population in V 2 largely exceeds the threshold 200 in the second and third year.

L =     0
We can also notice that a persistent colony appears in V 3 thanks to the individuals arriving from the other nodes.

Scenario 2: Time dependent distribution matrix

We consider a domain partitioned into three cells. The first one, V 1 , is located at a higher elevation than the two others, which are on the same level. We make the following choice of parameters : • In V 3 , we have a forest as in Section 4.1. Then we set α = 0.656, δ = 1, κ = 1.5 and take the same initial condition as in V 2 .

• In V 1 , α = δ = 1, κ = 1.
We assume now that during spring and summer, i.e. for t ∈ [0, 10], eventual dispersers tend to move to the second colony because of the abundance of food source in a crop landscape. While in autumn and winter, i.e. for t ∈ [START_REF] Johnson | Climatic warming disrupts recurrent alpine insect outbreaks[END_REF][START_REF] Thieme | Mathematics in Population Biology[END_REF],

they rather move downhill to the third cell, as we mentioned in Section 4.1. Therefore we consider the dispersal ratio matrix K 1 for t ∈ [0, 10] and K 2 for t ∈ [START_REF] Johnson | Climatic warming disrupts recurrent alpine insect outbreaks[END_REF][START_REF] Thieme | Mathematics in Population Biology[END_REF], In Figures 10 we present the evolution of the populations of voles in the three nodes for both models during three years. We can observe several dispersals occuring every year from each node. As we already noticed in previous simulations, the arrivals of voles increase the frequency of outbreaks. Moreover, we can see that in V 2 , that is the crop field landscape, the population of voles increases over the reproduction season with the appearance of peaks for the second model in Figure 10(d). However, due to the high mortality rate after harvesting, the increase of the total population each year is slow and contained, in particular its minimal value in a one year period is almost constant. On the contrary, in V 3 , the minimal value of the total population increases over one year, which means that there are more and more adults. At the end of the third year, we have approximatively 110 and 96 adult individuals in Figure 10(e) and Figure 10(f ) respectively. As a consequence of the increment of the adult population, we can see in the two figures that there are more and more dispersal outbreaks. Figure 10(f) also shows more peaks than in Figure 10(d).

K 1 =     0

  0165 females per female and per day during the reproduction season and only juvenile and adult individuals reproduce. In the time scale of our models this correspond to 0.3 per juvenile/adult vole per time unit. In Section 3.1, to validate our numerical schemes, we use as reproduction rate the piecewise constant function b(a) = 0.3χ [A 1 ,A 3 ](a) . Everywhere else in the paper we consider a smooth approximation of b, namely b

  The reproduction rate b 1 for α = 1. The mortality rate d 1 for δ = 1.

Figure 1 :

 1 Figure 1: The reproduction and mortality rates as functions of age.

  GD model, persistence.

  ID model, persistence.

  Both models, extinction.

Figure 2 :

 2 Figure 2: With reference to Section 2.3: the evolution of the vole population at one node. Left and center: low mortality in the GD and ID models respectively. Right : high mortality.

(3. 29 )

 29 Solution at V 2 before the arrival of voles from V 1 , t ∈ [0, ]. The dynamics is the same in the two models. The explicit form of the solution is given in(3.24). Computed solution for GD model for node V 2

Figure 3 :

 3 Figure 3: With reference to Section 3.1.1: the comparison between the population of voles in the exact solution and the GD model for η = 0.25.

Figure 4 :

 4 Figure 4: With reference to Section 3.1.1: the comparison between the population of voles in the exact solution and in the ID model.

2 Figure 5 : 2 Figure 6 :

 2526 Figure 5: L 1 -discrete relative errors for the GD Model in log/log scale at T = 0.6 and T = 1.2

  ),(3.18) and we compute the evolution of the vole populations in the two nodes by the GD model when the parameter of minimal time span for a dispersal is reduced to η = 3 × 10 -3 , c = 1 and the final time is T = 1.3. The result obtained for ∆a = ∆t = 3 × 10 -4 is in Figures 7.

Figure 7 :

 7 Figure 7: The voles populations of the two nodes computed with GD model for η = 3 × 10 -3

5 .

 5 Then we consider the travel time matrix and the distribution matrix as initial conditions are of 70 adult individuals in each node, equidistributed over their age class. The simulations, performed with ∆a = ∆t = 0.00625, are depicted in Figures 8.

Figure 8 :

 8 Figure 8: With reference to subsection 4.1: the evolution of the population over 4 years for the two models.

  5 and the initial condition consists of 100 adults individuals equidistributed over their age class. • In V 2 , we imagine to have a crop field as in Section 4.1. Then we set α = δ = 1, κ = 2.6 and the initial condition consists of 70 adults individuals equidistributed over their age class.

  the travel time between any two nodes is 0.5, so that

Figure 9 :

 9 Figure 9: With reference to Section 4.2.1: The population of voles over 3 years in V 1 , V 2 and V 3 for both models.

Figure 10 :Figure 11 :

 1011 Figure 10: With reference to Section 4.2.2: the population of voles over 3 years on three nodes representing different landscapes.

Figure 12 :

 12 Figure 12: Populations of voles in a colony starting with 19 juvenile (left) or 20 adult (right) individuals.

Figure 13 :

 13 Figure 13: Populations of voles in a colony starting with 18 juvenile (left) or 19 adult (right) individuals.

Table 1 :

 1 L 1 -discrete relative errors for the GD model

Acknowledgments

All the authors acknowledge the support of the Région Bourgogne Franche-Comté, projet 2017-2020 "Analyse mathématique et simulation numérique d'EDP issus de problèmes de contrôle et du trafic routier" and of the Zone Atelier Arc Jurassien projet 2018 "Modélisation et simulations numériques de la dynamique spatiale des populations de campagnols".

Conclusions

In this paper, we propose two models to describe the spatial dynamic of voles populations. The models are based on age-structured transport equations set on a graph. The nodes of the graph represent uniform landscape areas of relatively small size, where voles can settle. Edges are used to model the main feature of voles' macroscopic spatial dynamics, which is the dispersal of juvenile individuals whenever the total population at one node reaches a threshold number. This is supposed to be a gradual process in the first model while it takes place instantaneously in the second model. The parameters in the reproduction and mortality rates allows to reproduce the evolution in different landscapes.

The numerical simulations that we presented show the ability of the models to reproduce cycles of dispersals in the populations of voles, as well as persistence or extinction of a colony. The results that we obtained in this paper allow us to pursue the work with the construction of a hybrid model combining the ID model to describe the voles population dynamics at the nodes and the agents-based model developed in [START_REF] Marilleau | Coupling agent-based with equation-based models to study spatially explicit megapopulation dynamics[END_REF] to reproduce the spatial dynamics on large graphs. In a different direction it is also possible to enhance the two models by considering a distribution matrix which is not given a priori, but depends on the current state of the system. For example, one could imagine that dispersers would avoid an already overcrowded node.

Appendices A Allee's effect : survival number of a colony

We investigate in this part the minimal number of juvenile or adult individuals needed to maintain the existence of a colony when we consider the reproduction and the mortality rates defined by (2.9), (2.10) and (2.12) with α = δ = 1 and κ = 1.5. To that end, in the boundary condition at a = 0, (2.6), we consider a non constant function A of the form

so that the model can reproduce the Allee's effect. For simplicity, we just consider case in which β = 1, γ ∈ N, and

where φ * is a given value. As an example, if we take γ = 8 and φ * = 20 we obtain