Hafid Kadi

Mohammed Rebbah

Boudjelal Meftah

Olivier Lézoray

A Data Representation Model for Personalized Medicine

Keywords: Electronic health records, Personalized medicine, Data representation, Time-series data, Medical event, Clustering

Personalized medicine generates and exploits the patient data, such as genetic compositions, key biomarkers, treatment history, environmental factors, behavioral preferences, and demographics data. The information loss in the transformation process, the data types heterogeneity and the events time series in this set pose a problem in their exploration process. To solve these problems, we propose a data representation model. It considers the structured, temporal and/or non-temporal data and their types "Numeric, Nominal, Date, Boolean". After the data types "Date and Boolean" transformation, we treat the nominal data by dispersion and we apply several clustering techniques with different clusters numbers to control the numeric data distribution. Our work results in three homogeneous representations, these representations have only two dimensions and are easy to explore. Compared to the Symbolic Aggregate Approximation (SAX) technique, our model preserves time-series information, keeps data as much as possible and offers multiple simple representations to explore.

INTRODUCTION

Personalized medicine (PM) refers to the individualization of medical treatments based on the unique dataset of each patient. It generates and exploits stored patient data, which are often captured digitally in an "Electronic Health Record (EHR)" comprising profiles of many different patients. Essentially, an EHR refers to a longstanding, comprehensive health database resource that stores and manages all patient data files digitally under the custody of a licensed health entity. More specifically, it provides a digitalized view of the patient's demographics, data associated with the patient's clinical and medication history, diagnostic trajectory, social and economic environmental conditioning, geographical relocation, if any, as well as the patient's genetic data, if these exist [START_REF] Jensen | Mining electronic health records: towards better research applications and clinical care[END_REF].

Together, this massive data resource available via the EHR often includes not only homogeneous, heterogeneous, structured, unstructured and/or semi-structured data, but also the temporal and nontemporal data. Mixed in this huge bag of patient data are many captured medical events of different individual patients such as their body temperature measurements, blood pressure recordings and other time-series information, with different sorts and forms of data. As [START_REF] Ghazi | Mitigating the influence of the curse of dimensionality on time series similarity measures[END_REF] noted, we consider timeseries data to include all the observational sequences of a patient being captured vis-à-vis a medical event. Moreover, the EHR data resource contains a lot of hidden information and knowledge waiting to be mined and/or discovered. The process of reporting, evaluation, and medical decision-making based on the EHR data involves the extraction of relevant information and knowledge via specialized methods known as data mining techniques. The quality of information processing and knowledge discovery are thus directly linked to the availability, accessibility, type and form of the data to be extracted and aggregated for analysis. The objective of our work is to produce a high fidelity model for the representation of PM structured data. This is a challenging problem and our proposed model addresses several important scientific gaps: data heterogeneity, loss of data during data transformation, and interpretability of the representation over the course of a data mining process. To accomplish this non-trivial task, we represent the data by two parts. The first is dedicated to the representation of numeric data with clustering techniques, whereas the second part considers the representation of nominal data with respect to its dispersion. These two bodies of information are then joined into a single global representation table. Thanks to the simplicity of the obtained representation, healthcare specialists will be able to identify in the dataset both the key patient events, as well as the variations in the information conveyed by the data series. However, it is intended for the obtained representation to be used within automated medical decision-making processes such as disease prevention and/or adverse drug events prediction. Importantly, this paper emphasizes the need to explore the EHR data mining process that informs and challenges PM, which will ultimately enhance the ability of physicians and other care professionals to personalize high quality care to the inflicted individuals.

The rest of the paper is organized as follows. Section 2 explains the time-series representation process limitations. A novel data representation model proposed for PM is then detailed in Section 3 with Section 4 continuing on the discussion about the experimentation and the evaluation of the proposed model and the results analysis. The final section, Section 5, will provide concluding remarks and offer insights into practical implications and potential future works.

EHR INFORMATION & DATA LOSS

Most EHR data exploration models require a lot of the stored data, including temporal data, to be represented and transformed into an appropriate, meaningful and interoperable form. Research along these lines constitutes a common point among several past efforts. [START_REF] Bagattini | A classification framework for exploiting sparse multi-variate temporal features with application to adverse drug event detection in medical records[END_REF] proposed an approach that belongs to one of the medical branches advocating PM: the prediction of adverse drug events. The authors have used three phases: symbolic data representation, subsequences generation and classification. Their approach exploits sparse time series features, that is, it processes only the numeric data and ignores the other types of data that may appear in the EHR data. In the representation phase, they applied SAX technique in order to produce symbolic representation sequences for all the time-series. Regardless of the limits of the SAX technique and of the processed data qualities (different types), this approach provides a timely example on the use of PM data and puts forward the need for new techniques for EHR data representation. Among the emerging data mining methodologies, for example, the Deep Integrated Prediction (DIP) project [START_REF] Milad | A Predictive Approach Using Deep Feature Learning for Electronic Medical Records: A Comparative Study[END_REF] uses deep learning to represent numeric type attributes and the GloVe algorithm for discriminating nominal attributes in EHRs of patient records to showcase the cardiac disease. Graph-based Attention Model (GRAM) is an approach based on marking each medical visit events on a directed acyclic graph [START_REF] Edward | GRAM: graph-based attention model for healthcare representation learning[END_REF]. The work of [START_REF] Mallick | Entropy-based fuzzy hybrid framework for gene prediction network -an application to identify and rank the biomarkers for human lung adenocarcinoma[END_REF] presents an approach to study the genes interdependence in cancer cases. These researchers apply fuzzy logic to compute the information gain to represent the genes data vis-à-vis a graph. In [START_REF] Anima | Incorporating temporal EHR data in predictive models for risk stratification of renal function deterioration[END_REF], windows fixed with marking according to the time factor during diagnosis ('time-Windows') have been proposed. The time period was initially fixed to six months, thereby causing possible information loss concerning changes in the event's behavior. Years later, [START_REF] Jing | Learning from heterogeneous temporal data in electronic health records[END_REF] argue for a different approach to capture data temporality treatment as applied to the clinical events and representation extraction. For the transformation process, they used the Symbolic Aggregate Approximation (SAX) method on all observations to generate a representation as character chains, a process that may also cause the loss of information.

Treatment unification on the numerical data of patient records sometimes implies the need for data normalization. From a data transformation perspective, normalization assists in unifying the treatment scales so as to make the results comparable. Yet, this process can also generate information loss as encapsulated in the data series. Representation techniques applied on the results of these transformations can therefore lead to loss in the series as constraints such as the minimum length of the series to be processed may not be accommodated. In best practices, especially with time-series data, only one or two data types may generally be treated at the same time with the information and data loss problem omitted. To overcome challenges faced with the representation problems of time-series, data types homogeneity, and the information loss, an alternative data representation model is proposed. The proposed data model can represent the structured, temporal and/or non-temporal data and their differing types, including numeric, nominal, date, and Boolean. The temporal observations exist in a three-dimensional form (3-D: Patients, Events, Time). As an endeavor to unify the treatment process for the non-temporal data that are actually found in two-dimensions (2D: Patient, Data), we consider them in a three-axis form (Patient, Data, Time-1) such that "Time-1" takes only one value, i.e., "1". Such a model conserves as much data as is possible even for short time-series. Applying several techniques with a different number of clusters reallocate the data distribution. Accordingly, this work will result in a simplified representation with only two dimensions, making it easy to explore. With such results, it is clear that attention has now been given to the annotated challenges, including data type heterogeneity, the maximum coverage of data resource during treatment, the minimization of loss information and data during the transformation process.

For large-scale data transformation, imagine having a growing number of multiple data types to be extracted simultaneously. This would surely weigh down on the data exploration process, which is needed for the investigator(s) to obtain useful insights from the entire data set; more precisely, with satisfying the objective of having the emergence of irrelevant data to be minimized. Sometimes, novel data exploration approaches are used to represent data and extract dataset features, and may be enhanced with a combination of classification techniques and optimization tools such as the artificial bee colony (ABC) algorithm [START_REF] Razmjooy | A new design for PID controller by considering the operating points changes in Hydro-Turbine Connected to the equivalent network by using Invasive Weed Optimization (IWO) Algorithm[END_REF][START_REF] Khalilpour | Designing a Robust and Adaptive PID Controller for Gas Turbine Connected to the Generator[END_REF], an optimization technique based on variance reduction of the Gaussian distribution [START_REF] Namadchian | A New Meta-Heuristic Algorithm for Optimization based on Variance Reduction of Gaussian Distribution[END_REF], an invasive weed optimization approach [START_REF] Razmjooy | An improved quantum evolutionary algorithm based on invasive weed optimization[END_REF], an imperialist competitive algorithm [START_REF] Razmjooy | A hybrid neural network Imperialist Competitive Algorithm for skin color segmentation[END_REF], a particle swarm optimization (PSO) methodology [START_REF] Kolekar | Clinical decision making using artificial neural network with particle swarm optimization algorithm[END_REF] and/or a genetic algorithm (GA) [START_REF] Nguyen | Classification of Healthcare Data using Genetic Fuzzy Logic System and Wavelets[END_REF]) so as to train a classification and decision-making model. [START_REF] Yeh | A new hybrid approach for mining breast cancer pattern using discrete particle swarm optimization and statistical method[END_REF] and [START_REF] Nguyen | Classification of Healthcare Data using Genetic Fuzzy Logic System and Wavelets[END_REF], for example, are two medical data classification methodologies useful for EHR data exploration. The first separates patients with data foreshowing a high risk of breast cancer [START_REF] Yeh | A new hybrid approach for mining breast cancer pattern using discrete particle swarm optimization and statistical method[END_REF]. These researchers apply the PSO technique on features selection via statistical methods. The second uses the breast cancer and heart disease datasets for its exploration [START_REF] Nguyen | Classification of Healthcare Data using Genetic Fuzzy Logic System and Wavelets[END_REF]. The proposed model incorporates a feature extraction and transformation phase via the "Wavelet transformation" technique, followed by a training phase via the fuzzy standard additive model and the GA technique. The EHR data exploration can also involve the image-processing task stored in the data source. [START_REF] Kavya | Breast Cancer Lesion Detection From Cranial-Caudal View of Mammogram Images Using Statistical and Texture Features Extraction[END_REF] and [START_REF] Guo | A new interval differential equation for edge detection and determining breast cancer regions in mammography images[END_REF] are examples of medical image processing and breast cancer detection approaches. Symbolic Aggregate Approximation (SAX) is a representation and dimension reduction technique for time-series(e.g., [START_REF] Lin | Experiencing SAX: a novel symbolic representation of time series[END_REF][START_REF] Park | SAX-ARM: Deviant event pattern discovery from multivariate time series using symbolic aggregate approximation and association rule mining[END_REF]. It has the capacity to transform the numeric data series into sequences of consecutive symbols. Its first phase divides the time-series of length N into W equal segments where (W < N), and for each segment, the average of its values is computed. The second phase defines a corresponding strategy between a proposed alphabet with Z symbols and the new computed averages series. [START_REF] Jing | Learning from heterogeneous temporal data in electronic health records[END_REF] use the SAX technique to normalize time-series values. The normalization step generates two problems. To better understand these problems, we use two (2) time-series that represent the real temperature measures of two patients: X1={34.6, 34.6, 34.6, 34.6, 34.6, 34.0, 34.6, 34.6, 34.6, 34.6} and X2={36.6, 35.7, 35.5, 35.9, 36.1, 36.1, 35, 8, 36.8, 36.7, 36.0}. They will be considered to generate a string of nine characters at the base of three representation symbols. The first problem appears when applying the SAX method on X1 and X2 separately. Indeed, the generated representation is the same for both normalized v. non-normalized data. However the same symbol is used in the two representations that may correspond to two different intervals, for example, without normalization the technique SAX generates a symbol 'a'∈]-∞;34,462] for X1 and 'a'∈]-∞;35,939] for X2as shown in The second problem is observed after applying the multi-series SAX method, that is, processing over all time-series, based on the mean, variance and standard deviation or SD [START_REF] Razmjooy | A New Meta-Heuristic Optimization Algorithm Inspired by FIFA World Cup Competitions: Theory and Its Application in PID Designing for AVR System[END_REF]. The generated representations are not identical for normalized v. non-normalized series, but the behavioral observation of these series shows that there is a loss of meaning after the normalization process. Without normalization (Figure 2.a), the curve that presents the series X1 appeared totally under the curve of X2, in which all data values of X1 are inferior to X2 values. This apparently shows a meaningful behavior with the non-normalized dataset. Comparatively, with normalization (Figure 2.b), most, but not all, of the parties of the X2 curve are presented below the portions of X1, which implies the loss of data meaning and the behavioral direction of the two series. In addition to the information loss during normalization, the SAX technique requires an input parameter to specify the resulting symbol series length, which constitutes a third problem. This parameter must be less than or equal to the length of series to be represented; otherwise, all series that have a minimum length will be wasted. The impact of this criterion can lead to the loss of all data should and if the data contain only short series.

To eliminate these limitations, including the information and data loss, we propose a novel model representation that would be more protective on both information and data.

OUR PROPOSED MODEL

Problem Formulation

To formalize our proposed model, we use D to notify the EHR dataset that includes a patient set P={P1, P2, ..., Pn} and an event set E={E1, E2 , ..., Em} captured on time T as represented by the observational chronologies {T1, T2, ..., Tq} for each event.

Let eijr represents the observation value on the event Ei for the patient Pj and the chronology Tr. We then use Ei.T to present the longest series chronology in an event Ei. Let NuE represents the numeric events set and NoE represents the nominal events set:

where,

∀𝐸 ! ⊂ 𝑁𝑢𝐸˄∀𝑒 !"# ∈ 𝐸 ! ⇒ 𝑒 !"# ∈ ℝ
(1)

∀𝐸 ! ⊂ 𝑁𝑜𝐸˄∀𝑒 !"# ∈ 𝐸 ! ⇒ 𝑒 !"# is a String Type (2)
Equation 1 means that all observations of numerical events are real type values whereas Equation 2 means that all observations of nominal events are string type values.

We use IWki for intra-class inertia and IBki for inter-class inertia in relation to the event Ei and the clusters number ki, with the ITki to be their total inertia [START_REF] Choukri | Classification Algorithm Based on Nodes Similarity for MANETs[END_REF].

Description of the Model

Our approach treats structured, temporal and non-temporal data vis-à-vis their differing data types. As depicted in Figure 3,with data types to be inclusive of numeric, date, binary, and nominal, our Data Representation model per Region and Dispersion (DRRD) may be divided vertically into two main parts. The first part treats the numeric and date type data while the second part treats the nominal and Boolean type. Each part comprises several steps, including data representation and/or transformation, event marking, and result linearization. For all values ebj1 a transformation of the birth date event Eb is performed into the patient age. The other dates eijr are transformed by computing the observation appearance year according to the following formula:

Numeric & Date Data Representation

𝐷𝑎𝑡𝑒𝑇𝑜𝑁𝑢𝑚𝑒𝑟𝑖𝑐@𝑒 !"# 𝑒 $"% ⁄ B = 𝑁𝑏𝑌𝑒𝑎𝑟𝑠(𝑒 !"# -𝑒 $"% ) (3)
The function NbYears(eijr -ebj1) computes the difference between the date of eijr observation and the birth date ebj1 in terms of the years number, that is, it calculates the eijr observation age.

All simple or temporal data are presented as numeric temporal events, with this task to be presented as a time-series representation task. We consider event chronologies, and the result is a three-dimensional data reorganization (3D: P,NuE,T).

Partitioning Numeric Events

Clustering (or, Partitioning) is the task of grouping objects into subsets known as clusters based on a given similarity criteria between the object's properties to be examined [START_REF] Hancer | A new multi-objective differential evolution approach for simultaneous clustering and feature selection[END_REF]. Each clustering technique has a specific strategy for defining the clusters while associating an object to each cluster with respect to the cluster centers.

The clustering intends forming similar groups (clusters or classes) based on the ordered data of the same event. These clusters will be the belonging regions and the data comparison units of each patient on this event.

In this step (Part B1), for each event Ei in the reorganized data cube with time factor being neglected, we apply a clustering technique to produce ki clusters and their ordered centers list (C1, C2,..., Cp,..., Cu,...,Cki) as given in Equation 4. Figure 4details this step.

∀𝑝 < 𝑢˄∀𝑢 ≤ 𝑘 ! ⇒ 𝐶 & < 𝐶 ' (4)
The centers of clusters are ordered for their appropriate use with the symbols of an alphabet S to be defined later. Treatments sequentially include data selection, data sorting, data clustering, and finally the sorting of the resulting cluster centers as applied to each event Ei. The association of relevant data to the region of the nearest center will inform the data distribution in the regions on the ordered centers.

By distribution, we mean the manner of data dissemination, organization and display in the data presentation space. Four (4) clustering techniques have been evaluated: PAA or Piecewise Aggregate Approximation (e.g., [START_REF] Seunghye | Data representation for time series data mining: time domain approaches[END_REF][START_REF] Vineetha | An Analysis of Time Series Representation Methods Data Mining Applications Perspective[END_REF], K-MEANS (e.g., [START_REF] Hartigan | A K-Means Clustering Algorithm[END_REF], EM or Expectation-Maximization (e.g., [START_REF] Dempster | Maximum Likelihood from Incomplete Data via the EM Algorithm[END_REF], MDBC or Make Density Based Clusterer, based on hierarchical clustering techniques (e.g., [START_REF] Witten | Data Mining: Practical Machine Learning Tools and Techniques[END_REF]. Briefly, the principle of each of the referenced techniques in the sequel is as follows:

• PAA: Partitions a series of length n into N segments of length n/N, and the average data belonging to each segment is the representation of the latter.

• K-MEANS: Involves partitioning a dataset into a given number of clusters, each cluster being associated with a center point known as the centroid. Each point is assigned to the cluster with the closest centroid.

• EM: A probabilistic extension of the K-means algorithm, this iterative technique has been developed for incomplete data cases based on two steps; the first evaluates the expectation (E) while the second maximizes (M) the conceivable expectation.

• MDBC: A meta-clustering technique that relies on the results of a basic clustering approach to generate a probability distribution and a corresponding density for all observations. Specifically, the basic technique we use herein is the hierarchical clustering.

Each of these techniques will be applied to form ki clusters on each event Ei. Several ki values will be tested and the best result will be considered. Each event Ei must have different list values (Li) length with length(Li) greater than or equal to ki:

∀𝑒 !"# , 𝑒 !%( ∈ 𝐿 ! ⇒ 𝑒 !"# ≠ 𝑒 !%( ˄𝑙𝑒𝑛𝑔𝑡ℎ(𝐿 ! ) ≥ 𝑘 ! (5)
This is the task (Part C1) of notifying the values of each numeric reorganization cube event in a separate table. Each table generated for an event Ei has a number of columns, ColumnsNumber(Ei) being equal to:

𝐶𝑜𝑙𝑢𝑚𝑛𝑠𝑁𝑢𝑚𝑏𝑒𝑟(𝐸 ! ) = 𝑘 ! * 𝑙𝑒𝑛𝑔𝑡ℎ(𝐸 ! . 𝑇) (6)
where,

The function length(Ei. T) means the longest series T size in the event Ei.

The notification for each event Ei consists in indicating on its notification table MarkingTablei one cell for each value of each series. This indication must be pointed in front of the cell corresponding to the nearest cluster center, where,

∀𝑒 !"# ∈ 𝐸 ! , ∀𝐶 )* ∈ Z𝐶 + , … , 𝐶 , ! \ / 𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑒@𝐶 )* , 𝑒 !"# B = MIN &-+,…,, ! (𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑒@𝐶 & , 𝑒 !"# B) ⇒ 𝑀𝑎𝑟𝑘𝑖𝑛𝑔𝑇𝑎𝑏𝑙𝑒 ! [𝑗]e@(𝑟) * 𝑘 ! B + 𝑁𝐶g = X (7)
Distance (CNC,eijr) is the subtraction function between the CNC center and the eijr observation.

Numeric Data Marking

The marking process (see Algorithm 1) fills the cells indicated in the notifications table. We propose three (3) marking types. The first type is marking via real value as illustrated in Table 1, comprising the rewriting of the data value in the indicated cell (X=eijr). The second type, the binary marking (as illustrated in Table 2), is a presence signaling form. It notifies by "1" for the corresponding cells (X =1) and for the others by "0" (X=0). The last type is the marking via a symbol (as illustrated in Table 3). As centers are being ordered, we can associate a symbol from a predefined and ordered alphabet to each one. The marking task places the symbol that carries the same order with the center used during the notification (X = S[NC]).

Algorithm 1: Numeric Event Ei Marking.

The Numeric Data Linearization Task

Part D1 entails arranging the marking tables into a single table. Two (2) such arrangements are proposed:

• The first is the numeric event linearization (as illustrated in Table 4). In a table that collects all the events, we arrange(denoted by ∏) the marking tables one after the other,

where, • The second is the chronological linearization (as illustrated in Table 5). It joins (denoted by ∑ ) all columns of T1, the first times in a marking table with T2, the columns of the second times; then with T3 and so on until the last time is reached:

𝐸𝑣𝑒𝑛𝑡𝐿𝑖𝑛𝑒𝑎𝑟𝑖𝑧𝑎𝑡𝑖𝑜𝑛 = ∏ (𝑀𝑎𝑟𝑘𝑖𝑛𝑔𝑇𝑎𝑏𝑙𝑒 ! ) !-+,..,1 (8)
𝐶ℎ𝑟𝑜𝑛𝑜𝑙𝑜𝑔𝑖𝑐𝑎𝑙𝐿𝑖𝑛𝑒𝑎𝑟𝑖𝑧𝑎𝑡𝑖𝑜𝑛 = ∑ (𝑀𝑎𝑟𝑘𝑖𝑛𝑔𝑇𝑎𝑏𝑙𝑒 !2 " ) !-+,.

The Nominal & Boolean Data Representation

Part A2 considers the data of nominal and Boolean type.

The Boolean data transformation task converts temporal and non-temporal Boolean data into symbolic data. It entails replacing all values containing "0" by "F" and the value containing "1" by "Y":

m 𝐵𝑜𝑜𝑙𝑒𝑎𝑛𝑇𝑜𝑁𝑜𝑚𝑖𝑛𝑎𝑙@𝑒 !"# B = "F" if (𝑒 !"# = 0) 𝑎𝑛𝑑 𝐵𝑜𝑜𝑙𝑒𝑎𝑛𝑇𝑜𝑁𝑜𝑚𝑖𝑛𝑎𝑙(𝑒 !"# ) = "Y" if (𝑒 !"# = 1) (10)
The result will be a dataset that contains temporal as well as non-temporal nominal data. In this new nominal reorganization of events and following their chronologies, the data are represented according to three axes (P, NoE , T).

Dispersion "Scattering" of Nominal Events

Nominal Data Marking

The notification (Part C2) uses the result of the previous task. We use the chronology value r of each eijr observation as an index to notify all the nominal data on the dispersion tables:

∀𝑒 !"# ∈ 𝐸 ! , 𝐷𝑖𝑠𝑝𝑒𝑟𝑠𝑖𝑜𝑛𝑇𝑎𝑏𝑙𝑒 ! [𝑗][(𝑟 * 𝐿 ! . 𝑙𝑒𝑛𝑔𝑡ℎ) + 𝐿 ! . 𝑖𝑛𝑑𝑒𝑥𝑂𝑓(𝑒 !"# )] = X (12)
The function Li. indexOf(eijr) returns the eijr observation index in the set Li of different values to the event Ei.

The three (3) types of markings described above are reused. Algorithm 3 describes these steps.

The Nominal Data Linearization Task

Using the same numeric data linearization logic to collect the nominal events marking results into a single table (Part D2). There are always two (2) types of rearrangement, that via the nominal event v. that via the chronology of nominal observations.

Results Assembly

The two (2) linearization parts of numeric and nominal data require the assembly of their results. Two (2) proposals are advanced.

The first (part E1) joins the results of the numeric linearization and the nominal linearization in a single global representation for all patients. The consequences of this operation achieve a representation via real values, symbols, or binary representation. These results are two-dimensional tables, which highlights the simplicity of our final representation.

The second (part E2) is the view per patient (as illustrated in Table 6). Here, the proposal is to provide assistance to the health professionals, giving a clear idea of the variations in each event observations for a given patient. The view per patient comprises associating a data table to the patient Pj. The rows correspond to the captured events vis-à-vis the patient Pj representation. For each cell ki at each time Tr, we take only the non-null cells that contain the representation symbol. The created table must contain q columns so that Tq is the longest series Ei.T chronology for the event Ei. To provide the longest series Ei.T chronology information, we will color only the cells whose position is less or equal to the Ei.T length.

Table 6. View per patient prototype.

Events

T1 T2 T3 T4 ...

E1( On k1 Clusters ) E2( On k2 Clusters ) E3( On k3 Clusters )

... . ... ... ... ...

EXPERIMENTATION

For the model evaluation, we apply it to a real EHR dataset and compare the obtained representations with the results of the state-of-the-art SAX technique.

Dataset Description

We use the free data from OpenMRS Wiki (2018), system "Open Medical Record System". The OpenMRS platform in OpenMRS Project ( 2004) is an application for personalized EHRs. Based on 2,528 medical concepts, this dataset stores 476,973 observations for 5,284 patients. Two (2) diseases are observed in this dataset, namely, the "HIV Program" v. the "TB Program".

Data Selection/Transformation/Codification

For our experiments, we use only observations on the disease "TB Program". Only concepts that have more than one observation are selected. We obtained 21 numerical events, two (2) nominal events and four (4) date type events require a transformation operation. However, this dataset does not include observations of temporal nominal and Boolean type.

In the first step and after the transformation, the data have been re-arranged in a three-dimensional (3D) form, thereby repositioning twenty-five (25) numeric events.

Table 7 presents the statistics on ordered numerical events with each frequency. For the second step, the three-dimensional (3D) form of nominal events uses only the "GENDER" observations data and the "TRIBE" observations that are only available of this type.

Table 8 shows the statistics.

Results & Discussion

In line with the minimum length constraint (see Equation 5) vis-à-vis the list of different values applied on the numeric events to be partitioned, we have found that a single case generated by the event "980" has three (3) different values so that the clustering process is executed only for k= 2 and k= 3. The other cases of k = 4, 5 and 6 are not applicable.

The columns are named in relation to the event identifier, the cluster number and the observation number in the series, for example, if we take k =2 to partition the event "21 / HEMOGLOBIN," the maximum series of this event is observed on the patient number "3618" with three (3) values, which produces a representation with six (6) columns; in this case, "21_C0_S0, 21_C1_S0, 21_C0_S1, 21_C1_S1, 21_C0_S2, 21_C1_S2".

Table 9 summarizes the statistics of the lengths of series maximums and the produced columns number (i.e., the length of the empty notifications tables to be created) during the partitioning of each event. Data dispersion affects only the two (2) selected nominal concepts. The first "992843/GENDER" generates a dispersion table with two (2) columns"992843_0_M and 992843_0_F". The second concept "992844/TRIBE" generates a dispersion table with three (3) columns "992844_0_Luo, 992844_0_Luhya, 992844_0_Unknown".

Table 10 shows the data notification values for the numeric event "790/SERUM CREATININE" corresponding to the patient identified by "5126". We use the three (3) markings types, and the K-MEANS technique with k=2. The alphabet S={a,b} of marking symbols has only two (2) symbols vis-àvis the number of centers.

Table 10. Numeric event notification.

Also, Table 11 shows the nominal events "992843/GENDER" notification for the same patient. For each technique, we associate a best cases counter. This counter will be initialized to zero for each event; for each event and for each ki (ki=2,..., ki =6), we incrementally adjust the counter of the technique corresponding to the best inertia. If several techniques have the same best inertia in a given partition ki, we eliminate these results.

Altogether, Table 12 summarizes the global situation for the techniques used. For twenty-three (23) observed cases on twenty-five (25) events in the intra-class inertia evaluation, the K-MEANS technique is found to be superior to all the other approaches, followed by the MDBC, EM, and PAA techniques. The same result is obtained when evaluating the inter-class inertia, and for twenty-one (21) observed cases on twenty-five (25) events, the K-MEANS technique always remains the best.

Table 12. Global statistics of intra and inter class inertia.

The K-MEANS technique is superior among all others, and it can be used with all events. Yet, the qualities of each event, especially the distribution of the values, push us to choose the most appropriate classification technique for each event. To choose the technique and the number of clusters ki for each event Ei, we compute for all ki and all techniques the total inertia ITki:

𝐼𝑇 , ! = 𝐼𝑊 , ! + 𝐼𝐵 , ! (13)
The technique and the ki corresponding to the maximum ITki will be used for the event Ei representation.

Table 13 displays the evaluation results. For the nominal data dispersion statistics, no loss of information has occurred. The instances sum in each dispersion table of event "GENDER" (Figure 6.a) and event "TRIBE" (Figure 6.b) equals the total number of patients (824). Finally, the global representation has been generated on 431 columns. Of these, 426 columns associate with numeric events, several between them containing time-series in its initial form. The remaining 5 columns are associated with nominal events. These results prove the capacity of time-series treatments via our model.

Depending on the type of marking used, we obtain a global representation via real value (as shown in Figure 7), a binary global representation (as shown in Figure 8) or a global representation via symbol (as shown in Figure 9).

Representation Example & Evaluation

To evaluate our approach, we compute three (3) data symbolic representations of the patient identified with the id =75. Table 13 results have been enlisted for this purpose. The first representation demonstrates our DRRD approach whereby each event representation is taken separately from the others. The strategy entails concatenating the representation symbols of each event vis-à-vis their chronological order. The second and third representations apply the SAX technique with two (2) values different of the segments number W. W indicates the resulting representation length, which must be less than or equal to the series length in treatment. W =1 is used here as the minimum length acceptable by SAX in the second representation named SAX1, with W =10 as the maximum length corresponding to the longest series for the third representation named SAX10. The SAX technique requires an input symbol number Z comparable to the number of clusters. We then apply SAX1 and SAX10 on each event with Zi=ki.

Table 14 shows three (3) events symbolic representations to showcase the essential points in the comparison. Empty cells allude to the unobserved events (e.g., event 790) on this patient. Cells containing "/" are presented only for the SAX technique, and are informed by the minimum lengths conditions not respected for the series to be presented (e.g., event 21). For the nineteen (19) marked events, DRRD and the SAX1 generate new representations, whereas for the SAX10, new representations are generated for 8 events only. Additionally, our model represents each observation, whatever the length of the series. Even so, the SAX technique generates only the results with the same length if the series have passed the minimum length condition W. More generally, the SAX technique eliminates all series having a length less than W, which constitutes an information loss that directly affects the quality of the obtained representation. For instance, patients representation having only series lengths less than or equal to one will fail directly with the SAX technique and the parameter W≥2.

This example demonstrates the capacity of the proposed novel model to preserve time-series information and to conserve as much data as possible throughout the process of representation and transformation.

View per Patient Example

This solution consists in presenting a window on the data and its variations for each patient taken individually. Figure 10 shows an instance of this window for the patient identified by Id = 75. This data visualization can help practitioners in analyzing the observations. In this example, the patient has only one appearance of the "HEMOGLOBIN" event. The 'c' character encodes this event as it corresponds to the third level of this event among 4 levels (the number of clusters). The practitioner can also know the maximum length of the time-series of this event according to the number of colored cells. For the "HEMOGLOBIN" event, the time-series include a maximum of three samples. According to this information, a practitioner can easily decide if the patient has a strong need for other specimens or diagnostics. After certain experiments on the patients, the practitioners can notice then the key events for a given disease. This allows them to directly examine the desired event on this dashboard for the suspicious cases. Overall, this solution presents a detailed view of the data of a given patient in view of all the patient data.

CONCLUSION

In this paper, we first discuss PM, and its associated data, citing some types of patient profile data. The need to consider the composition of these data and to extract the hidden knowledge underlying the data have become a necessity for various critical data-intensive PM tasks, including reporting, evaluation, analysis, and medical decision-making.

Eventually, we argue that the patient data preparation and its representation can provide a staging basis for these different PM tasks. To date, the exploration of data for realizing PM has attracted a plethora of research. Notwithstanding, most of these works carry certain limits in their methodologies. Inevitably, the heterogeneity and number of data types, the information and data loss are the key challenges, and in order to showcase these difficulties, we cite some of the prominent works and describe certain limits in their results.

Subsequently, we propose a novel representation model working on numeric, nominal, date and Boolean data types to ensure maximum data resource coverage. By means of transformations and treatments on the nominal and numeric data, our proposed model overcomes the data types heterogeneity issues. The nominal data are processed by dispersion, and the numeric data is processed by data belonging region (Clustering and nearest centers). Such a resolution has a strong contribution on the time-series representation. This last is treated as a part of numeric data. Following the experiments, a single global representation is obtained for all patients. This representation includes the basic data detail even if the data have only one observation. It conserves the information embedded within the data especially for time-series as it represents each observation in relation to its original region.

On the individual patient scale, the results of the view per patient shows other details at the level of each patient, including the absence and presence of medical events, the number of tests and the specimens from the patient and their levels in relation to a global scale (number of clusters) of the event examined. The comparison of our approach vis-à-vis the SAX technique and its parameters shows the high efficiency and the capacity of our proposed model to minimize the information and time-series data loss (as compared to SAX10, where the segments number W is large) during the transformation and to provide strong descriptiveness of patient data (as compared to SAX1, where W is small). Put simply, our approach has the capability to process several data types at the same time whereas the SAX technique processes only the numeric data.

Our work simplifies the understanding of personalized medicine data by healthcare practitioners. It also provides a tool for the analysis and comparison of data between patients, and facilitates the tracking of key events and important observations by practitioners. Both diagnosis and treatment planning can benefit of such a tool. In particular, we believe that our proposed representation makes EHR data exploration easier and can help practitioners to turn towards personalized medicine as well as using computer-aided decision systems that mine EHR data.

The other essential point of our model is the three marking types proposal "By real value, Binary, By symbol", and the production of a simple, homogeneous and unified representation in the form of a table for each marking type. Achieving this solution will open the way for our future work, such as medical decision-making, classification and data exploration. Our representation richness in terms of data types, and data mining richness in terms of classification algorithms availability and their capacities for processing different data types, will give us the possibility to do other work that takes this representation as a basis for initial and forward learning. Currently, we are actively working to apply our representation model to other similar data resources, even if they are not medical. We believe the step ahead is to tackle the challenges in transforming and representing massive datasets containing time-series data.

  Figure 1.a. With normalization, it generates 'a'∈]-∞;0,771] for X1 and 'a'∈] -∞; 0,338] for X2as shown in Figure 1.b. The other symbols have the same behavior as the first symbol 'a'.

Figure 1 .

 1 Figure 1. SAX representation of: (a) non-normalized X1 and X2, (b) normalized X1 and X2.

Figure 2 .

 2 Figure 2. Multi-series SAX representation of: (a) non-normalized X1 and X2, (b) normalized X1 and X2.

Figure 3 .

 3 Figure 3. Proposed personalized medicine data representation model.

  This step selects the numerical type of data and transforms all date data types into a numeric (Part A1).

Figure 4 .

 4 Figure 4. An example of the clustering process of a numeric Event Ei.

  As argued by[START_REF] Chenguang | Development of a neuro-feedback game based on motor imagery EEG[END_REF], a better partition has either the lowest intra-class inertia, or the highest inter-class inertia. Hence, in order to compare the used classification techniques, we evaluate the results of our numerical events representation, at the base of intra and inter-class inertia. Radar graphs shown in Figure5.a and Figure 5.b have been generated based on the results of intra-class and inter-class inertia respectively. These graphs visually highlight the statistics of the used clustering techniques of PAA, K-MEANS, EM, and MDBC and illustrate the most appropriate clustering technique according to the event considered.

Figure 5 .

 5 Figure 5. Appropriate clustering techniques corresponding to statistics of: (a) intra-class inertia, (b) inter-class inertia.

Figure 6 .

 6 Figure 6. Dispersion of nominal events: (a) "GENDER", (b) "TRIBE".

Figure 7 .

 7 Figure 7. Part of a representation by real value.

Figure 8 .

 8 Figure 8. Binary representation part.

Figure 9 .

 9 Figure 9. Part of a representation by symbol.

Figure 10 .

 10 Figure 10. View per patient example "Patient Id=75".

Table 1 .

 1 Marking by Real value.

	EventMarking(P, Ei, Centersi, MarkingType){
	MarkingTablei=new Table[n][ki* Length(Ei. T)];
	if(MarkingType= Symbol)// Define an ordered alphabet S.
	Create an alphabet S of ordered symbols S0,S1,...,Sk;
	for all Pj in P do{
	int h=0;
	for all Tr in Ei.T do{
	NC=NearestCenter(Centersi , eijr); // Nearest center index
	if(MarkingType= Real Value)// Marking per real value.
	MarkingTablei[j][(h* ki)+NC]=eijr;
	if(MarkingType= Binary)	// Binary marking.
	MarkingTablei[j][(h* ki)+NC]= 1;
	if(MarkingType= Symbol)	// Marking per symbol.
	MarkingTablei[j][(h* ki)+NC]= SNC;
	h++;
	}
	}
	return MarkingTablei;
	}

Table 2 .

 2 Binary marking.

Table 3 .

 3 Marking by Symbol.

Table 4 .

 4 Event linearization.

	E1	E2	...	Em
	T1 T2 ... T1 T2 ... … T1 T2 …

Table 5 .

 5 Chronological linearization.

	.,1 "-+,…,3!4567(9 # .2)	(9)

.

  Nominal event dispersion process (Part B2) selects the different values set Li of each event Ei and creates a new corresponding empty table DispersionTableiWe associate Li. length cells for each value in the longest series Ei.T of this event Ei. The columns number in each table must be equal:∀𝑃 " ∈ 𝑃, 𝐷𝑖𝑠𝑝𝑒𝑟𝑠𝑖𝑜𝑛𝑇𝑎𝑏𝑙𝑒 ! e𝑃 " g. 𝑙𝑒𝑛𝑔𝑡ℎ = 𝐿 ! . 𝑙𝑒𝑛𝑔𝑡ℎ * 𝑙𝑒𝑛𝑔𝑡ℎ(𝐸

	}
	return DispersionTablei and Li;
	}	E1	E2 … Em	E1	E2 ... …	Em …
		T1	T1 … T1	T2	T2 … T2 …
	Algorithm 2 describes this step:
	Algorithm 2: Nominal Event Ei Dispersion.
	EventDispersion(Ei){
	Set Li = DistinctDataSelection(Ei); //Distinct data values selection.
	DispersionTablei=new Table[n][Li. length * length(Ei. T)];
	Y=0;
	For all z in (0,..., length(Ei. T)-1) do {
	For each value v in Li do {	//Labeling columns.
	DispersionTablei. Column[(z* Li. indexOf(v))+Y]. Name = Ei. ID+"_"+z+"_"+v;
	Y++;
	}

! . 𝑇) (11)

where, DispersionTablei

[Pj]

. length is the columns number in the table DispersionTablei with all patients Pj.

  Algorithm 3: Nominal Event Ei Marking.

	EventMarking(P, Ei, Li, DispersionTablei) {
	MarkingTablei= DispersionTablei;
	if(MarkingTypei= Symbol) {// Define an ordered alphabet S.
	w= Li. length-1;
	Create S of an ordered symbols S0,S1,...,Sw ;
	}
	for all Pj in P do {
	for all Tr in Ei.T do {
	SI= Li. indexOf(eijr);
	RI=(r* Li.length)+SI;
	if(MarkingType= Real Value)	// Marking per real value.
	MarkingTablei[j][RI]= eijr;
	if(MarkingType= Binary)	// Binary marking.
	MarkingTablei[j][RI]= 1;
	if(MarkingType= Symbol)	// Marking per symbol.
	MarkingTablei[j][RI]= SSI;
	}
	}
	return MarkingTablei;
	}

Table 7 .

 7 Statistics of numeric events.

	N°	Concept ID / Concept Name	Different instances number Frequency
	1	21 / HEMOGLOBIN	82	282
	2	654 / SERUM GLUTAMIC-PYRUVIC TRANSAMINASE	207	278
	3	678 / WHITE BLOOD CELLS	81	282
	4	729 / PLATELETS	204	282
	5	730 / CD4%	52	476
	6	790 / SERUM CREATININE	154	160
	7	851 / MEAN CORPUSCULAR VOLUME	58	282
	8	853 / CD8 COUNT	396	470
	9	952 / ABSOLUTE LYMPHOCYTE COUNT	190	275
	10 980 / BODY SURFACE AREA	3	3
	11 1113 / TUBERCULOSIS DRUG TREATMENT START DATE	7	9
	12 1279 / NUMBER OF WEEKS PREGNANT	19	38
	13 5085 / SYSTOLIC BLOOD PRESSURE	31	2. 144
	14 5086 / DIASTOLIC BLOOD PRESSURE	21	2. 143
	15 5087 / PULSE	113	2. 269
	16 5088 / TEMPERATURE (C)	64	2. 274
	17 5089 / WEIGHT (KG)	216	2. 262
	18 5090 / HEIGHT (CM)	79	135
	19 5092 / BLOOD OXYGEN SATURATION	25	2. 267
	20 5096 / RETURN VISIT DATE	2. 049	2. 265
	21 5242 / RESPIRATORY RATE	15	105
	22 5314 / HEAD CIRCUMFERENCE	31	122
	23 5497 / CD4 COUNT	330	478
	24 5599 / DATE OF CONFINEMENT	6	7
	25 5919 / BIRTH YEAR	803	824

Table 8 .

 8 Statistics of nominal events.

	N°	Concept ID / Concept Name	Different instances number Frequency
	1	992843 / GENDER	2	824
	2	992844 / TRIBE	3	824

Table 9 .

 9 Statistics of series maximums lengths and the columns produced number.

	N° Concept ID	Max Length of the longer series	Produced columns number k=2 k=3 k=4 k=5 k=6
	1 21	3	6	9	12	15	18
	2 654	2	4	6	8	10	12
	3 678	3	6	9	12	15	18
	4 729	3	6	9	12	15	18
	5 730	3	6	9	12	15	18
	6 790	2	4	6	8	10	12
	7 851	3	6	9	12	15	18
	8 853	3	6	9	12	15	18
	9 952	3	6	9	12	15	18
	10 980	1	2	3
	11 1113	3	6	9	12	15	18
	12 1279	4	8	12	16	20	24
	13 5085	10	20	30	40	50	60
	14 5086	10	20	30	40	50	60
	15 5087	10	20	30	40	50	60
	16 5088	10	20	30	40	50	60
	17 5089	10	20	30	40	50	60
	18 5090	4	8	12	16	20	24
	19 5092	10	20	30	40	50	60
	20 5096	9	18	27	36	45	54
	21 5242	4	8	12	16	20	24
	22 5314	4	8	12	16	20	24
	23 5497	3	6	9	12	15	18
	24 5599	2	4	6	8	10	12
	25 5919	1	2	3	4	5	6

Table 11 .

 11 Nominal event notification.

Table 13 .

 13 Chosen techniques and clusters numbers results.

N° Concepts ID Chosen technique Chosen ki

	21 5242		MDBC	2
	22 5314		K-MEANS	4
	23 5497		MDBC	2
	24 5599		PAA	3
	25 5919		MDBC	2
	Inertia	PAA	K-MEANS	EM	MDBC
	Intra class	4/23	8/23	5/23	6/23
	Inter class	0/21	10/21	3/21	8/21
	1 21		MDBC	4
	2 654		PAA	5
	3 678		EM	2
	4 729		MDBC	3
	5 730		EM	2
	6 790		PAA	2
	7 851		PAA	2
	8 853		MDBC	2
	9 952		EM	2
	10 980		PAA	2
	11 1113		PAA	4
	12 1279		MDBC	4
	13 5085		MDBC	4
	14 5086		MDBC	3
	15 5087		MDBC	4
	16 5088		MDBC	3
	17 5089		MDBC	4
	18 5090		MDBC	2
	19 5092		K-MEANS	6
	20 5096		MDBC	5

Table 14 .

 14 Representation by symbol of the patient identified by id = 75.

	Event Id Z=ki	DRRD	SAX1	SAX10
	21

ACKNOWLEDGMENTS

This work was realized as part of the Hubert Curien partnership (PHC) TASSILI cooperation program, between France and Algeria under the project code 19MDU212.


domain: info sdv

text:
Taha Boukhobza

Malika Smaïl-Tabbone

Automatic Synthesis of Boolean Networks from Biological Knowledge and Data

Keywords:

Boolean Networks (BNs) are a simple formalism used to study complex biological systems when the prediction of exact reaction times is not of interest. They play a key role to understand the dynamics of the studied systems and to predict their disruption in case of complex human diseases. BNs are generally built from experimental data and knowledge from the literature, either manually or with the aid of programs. The automatic synthesis of BNs is still a challenge for which several approaches have been proposed. In this paper, we propose ASKEeD-BN, a new approach based on Answer-Set Programming to synthesise BNs constrained in their structure and dynamics. By applying our method on several well-known biological systems, we provide empirical evidence that our approach can construct BNs in line with the provided constraints. We compare our approach with three existing methods (REVEAL, Best-Fit and caspo-TS) and show that our approach synthesises a small number of BNs which are covering a good proportion of the dynamical constraints, and that the variance of this coverage is low.

Introduction

Models of biological systems are important to understand the underlying processes in living organisms [START_REF] Lazebnik | Can a Biologist Fix a Radio?-Or, What I Learned While Studying Apoptosis[END_REF]. Once built, the model is an artefact that can be used to study a system through simulation. Several formalisms have been proposed to model biological systems [START_REF] Novère | Quantitative and Logic Modelling of Molecular and Gene Networks[END_REF], and they all have their own strengths and weaknesses. The choice of a formalism is guided by the question at hand: the best formalism is the most abstract formalism which can answer the question [START_REF] Bornholdt | Less Is More in Modeling Large Genetic Networks[END_REF]. For example, differential equations are a formalism suited to run detailed dynamic simulations because they contain information on kinetic parameters. However, they do not scale to large systems.

Boolean Networks (BNs) are a formalism used to study complex biological systems where prediction of exact reaction times is not of interest [START_REF] Albert | Boolean Modeling: A Logic-Based Dynamic Approach for Understanding Signaling and Regulatory Networks and for Making Useful Predictions[END_REF]. They play a key role to understand the dynamics of biological systems and predict their disruption in case of complex human diseases [START_REF] Biane | Abductive Network Action Inference for Targeted Therapy Discovery[END_REF]. The key notions of BNs are presented in Section 2.2. BNs are built from available knowledge about the structure of the system and data about the behaviour of its components (Section 2.3). The knowledge and data are used as constraints for the BN synthesis. The automatic synthesis of BNs from biological data and knowledge is still a challenge for which several methods have been developed. In Section 3, we review three state-of-the-art approaches: REVEAL, Best-Fit and caspo-TS.

In Section 4, we present ASKEeD-BN, a new automatic approach for the synthesis of BNs constrained in their structure and dynamics. We rely on the Answer-Set Programming framework to generate non-redundant BNs fulfilling the given constraints. We compare the performances of our approach with REVEAL, Best-Fit and caspo-TS on several biological systems with experimental and synthetic data (Section 5). Finally, we discuss the results and conclude.

2 Boolean Networks and their Synthesis

Prior Knowledge Network (PKN)

Part of the knowledge one has about a biological system is the list of components (genes, proteins. . . ) constituting the system and how these components influence each other. Influences have a polarity: activation (polarity "+") or inhibition (polarity "-"). The parents of a component X are the components which influence X. A Prior Knowledge Network (PKN) encodes this knowledge. The nodes of the network are the components of the system. The edges are directed from parent components to child components and labelled "+" or "-" according to the polarity of the influences. Fig. 1 shows an example PKN for a system of three components. In this PKN, C and A are the parents of C.

-"A activates C" -"B interacts with itself" -"C activates A" -"C interacts with B" -"C inhibits itself"

Boolean Networks (BNs)

BNs were introduced by Kauffman [START_REF] Kauffman | Metabolic Stability and Epigenesis in Randomly Constructed Genetic Nets[END_REF] to model genetic regulatory networks. Concepts used in BNs are described in a recent review [START_REF] Schwab | Concepts in Boolean Network Modeling: What Do They All Mean?[END_REF]. Two examples of BNs are given in Fig. 2.

The components of a BN are the components of the considered biological system. For example, a BN modelling a system of three proteins Like for the PKN, the structure of a BN is defined in terms of parent-child relationships between the components. A component P which appears in the transition function of a component X is called a parent of X. If the parent is negated in the DNF, we say that the polarity of the influence of P on X is negative. Conversely, if the parent is not negated, the polarity is positive. The Interaction Graph (IG) summarises these relationships as a directed graph. The directed edge P → X is labelled with "+" or "-" depending on the polarity of the influence P has on X. For example, the IG of B 1 contains B + -→ B and C --→ B because B appears positively and C appears negatively in the transition function associated to B. As we will see in Section 2.3, the PKN will act as a hard constraint on the IG of the BNs we want to synthesise.

The BN dynamics is obtained by applying iteratively the transition functions starting from each possible configuration. The order of application of the transition functions is defined by the update scheme. The synchronous, asynchronous and mixed update schemes are the most commonly used. In the synchronous update scheme, the transition functions are applied all at once, while in the asynchronous scheme, they are applied one by one. In the mixed update scheme, any number of components can be updated at each step. Thus, the update possibilities from both the synchronous and asynchronous update schemes are included in this third update scheme.

The State Transition Graph (STG) is a directed graph whose nodes are the 2 n possible configurations of the BN. In this graph, there is a directed edge from c to c if c is the result of applying to c the transition function(s) according to the chosen update scheme. Fig. 2 shows examples of synchronous and asynchronous STGs. As we will see later, the synthesis of BNs constrained in their dynamics implies to enforce their STG to contain specific edges, corresponding to specific transitions of configuration. We will also see how we use the mixed STG to quantify how the synthesised BNs match the dynamical constraints.

Synthesis of BNs from PKN and multivariate TS

In general, BNs that model biological systems have to satisfy two categories of constraints. On one hand, the BNs have to comply with a PKN. The PKN constrains the structure of the synthesised BNs by defining which components can appear as variables in each transition function and the polarity of those variables. Hence, a component P is allowed to appear in the transition function of a component X with a polarity s if the PKN contains an edge P s -→ X. Formally, a BN is compatible with a PKN if its IG is a spanning subgraph of the PKN. In other words, the IG of a BN compatible with a given PKN is formed from the vertices and a subset of the edges of the PKN. For example, the two BNs presented in Figs. 2a and2b are compatible with the PKN given in Fig. 1. On the contrary, a BN containing the transition function f A := B is not, since the IG of this BN contains the edge B + -→ A, which is not in the PKN. A BN having f A := ¬C is also incompatible: despite C being a possible parent of A, the negative polarity is not allowed, since the PKN does not contain the edge C --→ A. On the other hand, the synthesised BNs are expected to reproduce as well as possible the sequence of configurations extracted from an observed continuous multivariate Time Series (TS) of the concentration of the components over time. An example of multivariate TS is given in Table 1. Various strategies for extracting the sequence of configurations and fitting the transition functions to the observations were proposed in the literature.

We focus here on the automatic synthesis of BNs that respect the structure of a given PKN and are designed to reproduce as well as possible the observations from one given multivariate TS. For each synthesised BN, this ability of reproducing the observations is measured in terms of coverage proportion, i.e., the proportion of transitions observed in the multivariate TS that are retrieved by the BN when computing its STG according to the mixed update scheme. Ideally, an identification method would only return BNs with a perfect coverage proportion (i.e., 1).

Table 1. Multivariate TS of the three-components system given as example. The continuous concentrations of the components have been sampled for 20 time steps. Here, all the observations range from 0 to 100. The value resulting from the binarisation with a threshold of 50 is indicated by the colour of the cells: green if the result of the binarisation is 1 and red if 0. The resulting binary vectors are the configurations. Here there are four configurations (010, 011, 100 and 001) lasting Several studies have been dedicated to the automatic synthesis of BNs from PKNs and observed multivariate TS. Here, we review three main state-of-the-art approaches: REVEAL [START_REF] Liang | REVEAL, a General Reverse Engineering Algorithm for Inference of Genetic Network Architectures[END_REF], Best-Fit [START_REF] Lähdesmäki | On Learning Gene Regulatory Networks under the Boolean Network Model[END_REF] and caspo-TS [START_REF] Ostrowski | Boolean Network Identification from Perturbation Time Series Data Combining Dynamics Abstraction and Logic Programming[END_REF].

For each component of the system, REVEAL tests all the possible combinations of its parent nodes, and attempts to find the functions that explain all the observations of the binarised TS. For example with the multivariate TS from Table 1: REVEAL tries to explain 010 → 010 → 010 → 010 → 011 → 011 → 011 → 100 → . . . Hence, it cannot handle inconsistencies -such as a configuration being associated to distinct successor configurations. Such inconsistencies are frequent when sampling concentration along time because the processes implicated do not all have the same speed. In the example (Table 1), observing both 010 → 010 and 010 → 011 is an inconsistency which causes the failure of REVEAL. Furthermore, REVEAL cannot use the influence signs from the PKN, and since it uses an already binarised TS, it is possibly biased by the chosen binarisation.

Like REVEAL, Best-Fit tests every possible combination of the parent nodes of each component. It cannot use the influence signs and works on the binarised TS as well. Unlike REVEAL it can manage inconsistencies from the TS since it returns the functions that explain the maximal number of time steps. In Table 1, since 010 → 010 is observed three times and 010 → 011 only once, Best-Fit will focus on explaining the former instead of the latter.

caspo-TS was designed to manage several multivariate TS, corresponding to several experiments where the system is perturbed (forced activation or inhibition of some components), and where some measurements are potentially missing. Unlike REVEAL and Best-Fit, caspo-TS takes the influence signs into account. But it can only generate locally monotonous BNs, i.e., BNs for which a parent of a component cannot be both its activator and its inhibitor. B 2 is an example of BN caspo-TS cannot generate due to non local monotony. Indeed, in f B , the components B and C act both as activator and inhibitor of B. caspo-TS works as the following: first, it derives the set of BNs that are compatible with the given PKN and an over-approximation of the dynamics of the TS, using the so-called most-permissive semantics [START_REF] Chatain | Most Permissive Semantics of Boolean Networks[END_REF]. Because of this over-approximation, the result can contain many false positive BNs, i.e.,BNs optimising the cost function used under the hood of caspo-TS, while their asynchronous dynamics is not able to reproduce the configurations sequence of the multivariate TS. These false positive BNs are subsequently ruled out using exact model checking. This filtering is PSPACE-hard, but thanks to the first step, a large set of BNs has already been excluded. [START_REF] Chatain | Most Permissive Semantics of Boolean Networks[END_REF] Our Approach: ASKEeD-BN

Details of the Approach

We propose an approach for the Automatic Synthesis of Boolean Networks from Knowledge and Data (ASKEeD-BN). It computes a non-redundant set of BNs complying with a given PKN and one observed multivariate TS. Unlike REVEAL and Best-Fit, ASKEeD-BN is capable of using the influence signs provided in the given Prior Knowledge Network (PKN) and the raw values of the input multivariate Time-Series (TS). Unlike caspo-TS, ASKEeD-BN directly fits the behaviour of each component with the TS. Also, it is not limited to the synthesis of locally-monotonous BNs.

For each component of the studied system, our approach searches among all possible transition functions. All the transition functions that do not respect the given PKN are ruled out. Then, every remaining candidate is evaluated on the basis of both their simplicity and their ability to reproduce the given observations. The candidate transition functions for the component X might not be able to explain all the binary state transitions happening at time t → t . The set of unexplained t is denoted U. Every time step t in U is associated with a measure stating "how far" the continuous value X t is from the binarisation threshold θ: |θ -X t |. These spotted errors are then averaged on the T time steps of the TS through the Mean Absolute Error (MAE):

MAE X = t ∈U |θ -X t | T
Among the candidates having the smallest MAE, we select the ones that has the smallest number of influences. Finally, we create all the possible BNs by generating all the combinations of the selected functions. We implemented our approach using Python and the Answer-Set Programming framework (ASP) with the system clingo [START_REF] Gebser | Answer Set Solving in Practice[END_REF]. ASP is a declarative programming language oriented towards difficult (NP-hard) search problems. The possible solutions of a problem are described with the constraints they must fulfill. These constraints are written as a logic program. The ASP solver is tasked with finding the solutions of the program. To do so, it uses a Conflict-Driven Clause Learning (CDCL) algorithm inspired by SAT solvers. In our case, the CDCL algorithm avoids the evaluation of all the possible transition functions by learning from conflicts: whenever it finds that a candidate is in conflict with the constraints, it creates a new constraint that explains the conflict. These learned constraints subsequently eliminate other conflicting candidates, pruning the search space. Thanks to these pruning heuristics, our approach is efficient. ASP and in particular clingo, have already been used in similar contexts including caspo-TS.

Illustration on the Toy Example

Let us illustrate our approach on the toy example consisting of the PKN in Fig. 1 and the multivariate TS in Table 1.

When no PKN is available, the default PKN is a complete graph assuming that each component can inhibit / activate all the others (including itself). In this setup, a component with n parents have 2 2 n possible transition functions. In the toy example, each component can be explained by 2 2 3 = 256 distinct functions, which correspond to 16 777 216 potential BNs (formed by all the possible combinations of all the candidates of each component). Thanks to the available PKN, the number of candidate functions for each components A, B and C falls respectively to 3, 16 and 6. Besides the CDCL pruning, ASKEeD-BN virtually evaluates all the candidates, but for illustration purpose we will focus on the two that are present in B 1 and B 2 (Figs. 2a and2b).

For the component A, the candidate f A := 0 does not contain any literal and it cannot explain the transition of configuration for A at t 7 → t 8 . Hence, the set U of unexplained time steps is {t 8 }. The concentration of A at time t 8 is 61, and the candidate's MAE is thus |50-61|/20 = 0.55. The candidate f A := C involves one literal (which is C). This candidate can explain all transitions. Hence, U = ∅ and the MAE associated with this candidate is 0. Despite requiring more literals, f A := C is a better candidate than f A := 0 because its MAE is smaller. The comparisons of the candidates proposed for the components B and C in B 1 and B 2 are summarised in Table 2.

For the toy example, our approach returns B 1 as the only solution. It retrieves the 3 configuration transitions extracted from the binarised TS, thus its coverage proportion is 1. REVEAL does not find any BN, and the BN returned by Best-Fit does not comply with the PKN. caspo-TS finds 5 BNs with coverage proportions ranging from 0.33 to 1 (standard deviation of 0.25).

Table 2. Number of influences and MAE for the candidate functions in B1 (Fig. 2a) and B2 (Fig. 2b). A checkmark indicates the candidate selected by our approach, and the best for each criterion: (1) minimal MAE and (2) minimal number of influences. 5 Datasets and Procedure for the Comparative Evaluation

Datasets

In order to compare our approach with REVEAL, Best-Fit and caspo-TS, we used eight biological systems. For two of these systems (yeast's cell cycle and A. thaliana's circadian clock), their PKN and experimental multivariate TS are taken from [START_REF] Locke | Experimental Validation of a Predicted Feedback Loop in the Multi-Oscillator Clock of Arabidopsis Thaliana[END_REF] and [START_REF] Spellman | Comprehensive Identification of Cell Cycle-Regulated Genes of the Yeast Saccharomyces Cerevisiae by Microarray Hybridization[END_REF] respectively. These two systems are summarised in Table 3. They respectively involve 4 and 5 components. For the six other systems 3 , we conducted our experiments on multivariate TS that we simulated from existing BNs taken from the repository of example BNs of the package PyBoolNet [START_REF] Klarner | PyBoolNet: A Python Package for the Generation, Analysis and Visualization of Boolean Networks[END_REF]. For these systems, the number of components ranges from 3 to 10. For each system, the used PKN is the IG of the associated BN. As for the generation of the multivariate TS, three parameters are taken into consideration: the update scheme (in {synchronous, asynchronous}), the maximum number of introduced repetitions of each configuration (in {1, 4}) and the standard deviation of the added noise (in {0, 0.1}). For each setting of these parameters, we follow a procedure similar to what is implemented in the generateTimeSeries function of the R package BoolNet [START_REF] Müssel | BoolNet-an R Package for Generation, Reconstruction and Analysis of Boolean Networks[END_REF]:

1. choose randomly a configuration of the considered BN, 3 raf, randomnet n7k3, xiao wnt5a, arellano rootstem, davidich yeast and faure cellcycle 2. on this configuration, apply the update function(s) 20 times w.r.t the chosen update scheme, 3. duplicate randomly each configuration in the obtained sequence (added in contrast to generateTimeSeries), 4. add a Gaussian noise with a standard deviation of N.

For a given setting of the 3 parameters and a given system, we run the procedure 7 times (with different random seed). In the following, we denote ARN the setting with the Asynchronous update scheme, Repetitions (of 4) and Noise (of 0.1). We believe that this setting allows us to obtain multivariate TS which are quite close to real TS.

We illustrate here how to generate a synthetic multivariate TS in the ARN setting for B 1 (Fig. 2a). We would start from a random configuration. Let it be 010. Then we apply 20 times the transition functions of B 1 with the asynchronous update scheme. This process is not deterministic as any path from Fig. 2e (right) starting from 010 and of length 20 is valid. Let's say we obtain a path starting with 010 → 011 → 010 → 011 → 111 → 101 → . . . Then we add a random number of duplications (in bold). The beginning of the sequence could for example look like 010 → 011 → 011→ 010 → 011 → 011→ 011→ 111 → 101 → 101→ 101→ 101→ . . . Finally, we add a random Gaussian noise with a standard deviation of 0.1. The synthetic multivariate TS could now start with (0.02; 0.92; -0.16) → (0.04; 0.8; 0.7) → (-0.05; 1.06; 0.7) → . . .

Details on the Evaluation Procedure

For REVEAL and Best-Fit we use the implementation from the R package BoolNet [START_REF] Müssel | BoolNet-an R Package for Generation, Reconstruction and Analysis of Boolean Networks[END_REF]. caspo-TS is ran with the option mincard, that asks for BNs with functions minimising the number of influences. Note that this is also what our method optimises.

In the following, we define an experiment as a BN identification method applied on a system with one multivariate TS. The unicity of the multivariate TS makes the problem under-specified and allows us to evaluate the performances of the different approaches in this context.

REVEAL, Best-Fit and our approach need the binarised multivariate TS in their inputs. We use a simple form of binarisation: the binarisation threshold is defined as min + (max -min)/2. All values from the multivariate TS greater or equal to the threshold are binarised to 1, and to 0 otherwise. For the two systems with real TS, the theoretic range of the values is not know in advance, so the binarisation threshold is determined component-wise: the components are binarised taking into account their observed minimum and maximum. For the six systems with the synthetic multivariate TS, we know a priori that the values of all the components are between 0 and 1 (± the noise). In case of noisy data, the fluctuations of a constant component are interpreted as state change when using a threshold computed component-wise. However, the identification methods are not capable to detect these spurious transition in the binarized TS.

Hence, we compute the binarisation threshold globally, on all the observation of all the components.

In order to have a fair comparison of the methods, and since caspo-TS is making the binarisation itself and is not aware that the theoretical minimum and maximum of the components are 0 and 1 (± the noise), we correct a posteriori the transition functions it returned. The value of the constant is set to the binarised value that is the most present in the binarised TS of the component concerned. Also, since caspo-TS does not return a function for the components without parents in the PKN nor for the components that it founds constant for all the TS (in the case where no noise is involved), we use the same technique to set the transition functions to their correct values. We also added a step to filter out BNs returned by REVEAL and Best-Fit which do not respect the polarities given in the PKN.

For all the BNs returned by the four methods (and after the PKN-based filtering for REVEAL and Best-Fit), we use PyBoolNet [START_REF] Klarner | PyBoolNet: A Python Package for the Generation, Analysis and Visualization of Boolean Networks[END_REF] to compute the STG of each retrieved BN according to the mixed update scheme. Finally, we evaluate the results of each experience according to three criteria:

the number of BNs returned; the median of the coverage ratios: the proportion of configuration transitions extracted from the input TS that are present in the mixed STG; the standard deviation of the coverage ratios.

All data and programs needed to reproduce the presented results are accessible at https://gitlab.inria.fr/avaginay/OLA2021.

Results

Results on Systems with Real PKN and Experimental

Multivariate TS Yeast (Fig. 3 left). For this system caspo-TS find 61 BNs while Best-Fit and ASKEeD-BN both find 16 BNs. As for REVEAL, due to inconsistencies in the TS, it does not return any BN. Concerning the coverage, on the 7 transitions observed in the TS, the BNs synthesised by Best-Fit recover 4 and the BNs synthesised by ASKEeD-BN recover five. The best coverage ratio (6 retrieved transitions over 7) is obtained for 8 BNs synthesised by caspo-TS (among the total of 61). Nevertheless, as the box plot shows, the BNs synthesised by caspo-TS present a large variance in their coverage.

A. thaliana (Fig. 3 right). For this system, REVEAL returns no BN. The only BN returned by Best-Fit has all the components set to 1 and recovers 5 transitions over the 10 observed. ASKEeD-BN also returns a single BN with a perfect coverage since the BN recovers all the 10 transitions. As for the 5 BNs synthesised by caspo-TS, we can make the same observation as before: they present a variability in their coverage. The best coverage obtained by caspo-TS are from 2 different BNs including the one synthesised by ASKEeD-BN. To sum up, the results on these two real examples show that: -REVEAL constantly fails to return any BN. At the opposite, caspo-TS returns more BNs than the other methods; the coverage of the BNs returned by both our approach and caspo-TS are better than for Best-Fit; caspo-TS presents worse variability in the coverage ratio of its BNs compared to our approach.

Results on Systems with generated multivariate TS

Number of Synthesised BNs: The total number of BNs returned on the synthetic datasets and the number of times the identification methods failed returning any BNs are reported in Table 4. The table shows that a large proportion of BNs generated by REVEAL and Best-Fit were not complying with the influence signs from the input PKN. The following reported results do not take into account these non-compliant BNs. REVEAL is the method which returns the smallest number of BNs, in particular in the ARN setting. This is due to the inconsistencies in the TS, which are frequent in the ARN setting (as in real TS). On the opposite, caspo-TS is the method that returned the largest number of BNs. Moreover, when considering all experiments, there are 18 experiments for which caspo-TS generated more than 100 BNs. In these cases, we stopped the enumeration and analysed the 100 first BNs caspo-TS returned. Despite this limit, caspo-TS returned between 5 and 7 times more BNs than our method.

From here on, we focus on the results of the experiments corresponding to the ARN setting (Asynchronous update scheme, random Repetition of configurations, and Noise addition) after having remove the BNs from REVEAL and Best-Fit which does not respect the given PKN.

Coverage ratio: To assess the coverage ratio criterion, instead of plotting the boxplots for the 42 experiments of this setting (6 systems times 7 replicates), we summarised them in Fig. 4. In the scatter plot, each experiment is represented by a point whose coordinates are the coverage ratio median of the synthesised BNs and the associated standard deviation (std). The more top-right a point is, the better the corresponding identification method is (i.e., it produces BNs with high coverage ratio and low std). We can see that for the few experiments for which REVEAL was able to return BNs, the median coverage is actually excellent. The median coverage of the BNs returned by Best-Fit is almost uniform: Best-Fit lacks regularity in finding BNs with good coverage. But the high pick around 0 on the plot of std distribution shows that for a given experiment, the BNs returned by Best-Fit have similar coverage rates. caspo-TS and our approach have a very similar distribution of median coverage. They are both good at finding BNs with very good coverage. But here again, for a given experiment, the BNs synthesised by caspo-TS present a bigger variation of their coverage proportions than the ones synthesised by our approach.

Conclusion and Perspectives

We presented ASKEeD-BN, a novel method to create BNs from a PKN and a multivariate TS. The results on 8 biological systems showed that our approach has the best trade-off on the evaluation criteria: it returns a small set of BNs with a high coverage median and low variance. Our results actually confirm that although caspo-TS finds good BNs, too many sub-optimal BNs are also retrieved. Indeed a new version of caspo-TS was recently proposed to tackle this problem [START_REF] Chevalier | Synthesis and Simulation of Ensembles of Boolean Networks for Cell Fate Decision[END_REF].

We now present two perspectives to improve our approach and the study. First of all, real datasets may contain outlier measurements which could mislead the computation procedure of the binarisation thresholds we used in this paper. It would be interesting to see how such cases impact the performances of the identification methods and to propose a better binarisation procedure with prior outliers detection for instance. Second, contrarily to REVEAL, Best-Fit and caspo-TS, our approach does not handle multiple multivariate TS. However, biologists often have several multivariate TS generated with perturbations forcing some components to stay either active or inactive. However, exploiting such Fig. 4. On the scatter plot, each point represents an experiment in the ARN condition for which a given method potentially returned several BNs with different coverage ratios. The horizontal coordinate of the point is the median of these ratios. The vertical coordinate is their standard deviation (std). For a better visualisation, the coordinates have been jittered with a variance of 0.1 on both axes. The curves on the top (resp. on the left) of the scatter plot are the probability densities of the median (resp. the std) of the points in the scatter plots. The densities have been estimated from the non-jittered coordinates of the points with the Gaussian kernel density estimation method. The smoothing parameter of the estimator was determined automatically (with the Scott method). The areas under all these curves are 1, and the picks show where the points are the most concentrated.

supplementary data gives more information about the behaviour of the studied system in specific conditions (e.g., pathological states). This knowledge allows to constrain even more the space of solutions.

Finally, we are currently working on an automatic pipeline for BN synthesis from a curated mathematical model repository, namely BioModels [START_REF] Malik-Sheriff | BioModels-15 Years of Sharing Computational Models in Life Science[END_REF]. This requires (i) automatic extraction of the PKN from the model structure encoded in the SBML4 file format and (ii) generation of a multivariate TS by simulation of these models.

Fig. 1 .

 1 Fig.1. PKN example of a three-components system.

Fig. 2 .

 2 Fig. 2. The transition functions, derived interaction graph, and state transition graphs according to synchronous and asynchronous update schemes of two BNs.

  is a vector which associates a Boolean value (1/active or 0/inactive) to each component of the BN. A BN with n components has 2 n possible configurations. For example, the 2 3 = 8 possible configurations of a BN with 3 components are: 000, 001, 010, 011, 100, 101, 110 and 111. Each component has an associated transition function (B n → B) which maps the configurations of the BN to the next value of the component. The transition functions are usually written as Boolean expressions. In this paper, these expressions are in Disjunctive Normal Form (DNF), i.e., disjunctions of conjunctions. The conjunctions are satisfiable, which means they do not contain a literal and its contrary. The operators ¬, ∧, ∨ represent respectively negation, conjunction and disjunction. Figs. 2a and 2b show examples of transition functions. The transition function associated to B in B 2 states that the value of B will be 1 if either the value of B or of C was 1 in the previous configuration.

3

  respectively 4, 3, 3 and 10 time steps. Vertical bars indicate a change of configuration. configurations sequence: 010 → 011 → 100 → 001 time 1 State-of-the-Art Methods of BN synthesis from PKN and TS

  candidate fB := B ∧ ¬C fB := (B ∧ ¬C) ∨ (¬B ∧ C) fC := ¬C fC :

Fig. 3 .

 3 Fig. 3. Number of transitions retrieved by the BNs synthesised using the different methods on the systems yeast (left) and A. thaliana (right). The blue dashed line indicates the number of transitions that were observed in the multivariate TS.

Table 3 .

 3 Summary of two biological systems and their corresponding datasets

	System	Genes				PKN				TS	Source
	yeast (cell cycle)	Fkh2, Swi5, Sic1 & Clb1	Sic1 does not influence itself nor Fkh2	14 time steps 6 transitions	[18]
	A. thaliana (circadian clock)	LHY, PRR7, TOC1, X & Y	X	+	TOC1 + -	LHY +	-	-	Y -+	PRR7	50 time steps 11 transitions	[13]

Table 4 .

 4 Number of experiments for which each method failed to return any BN, number of BNs returned over all 336 experiments with synthetic TS and number of BNs returned over the 42 experiments with the ARN setting. The "before / after filter" refers to the step which rules out the BNs not respecting the signs of the given PKN (see Section 5.2).

	REVEAL setting before filter after filter before filter after filter Best-Fit	caspo-TS ASKEeD-BN
	# failing experiments all	230	240	0	64	20	0
	# total returned BNs all 100 677 500	406	100 678 198	724	8481	1210
	# total returned BNs ARN	3	3	51	35	720	85

Systems Biology Markup Language

Acknowledgements We thank Julie Lao and Hans-Jörg Schurr for their valuable comments and suggestions.