Émilien Gauthier

Laurent Brisson

Philippe Lenca

Franc ¸oise Clavel-Chapelon

Stéphane Ragusa

Challenges to build an information system for a breast cancer risk score

Cancer has become the leading cause of death worldwide according to the World Health Organization. As a consequence, prevention and screening programs are set up by health authorities to decrease its incidence. Programs efficiency can be increased by targeting highest risk subset of the population. Efficient information systems capable of monitoring the population risk are thus needed. Constraints to build such a cancer risk score and their impacts on the information system are presented. As it will be shown, beyond risk score performance, a major constraint concerns the place of domain expert and the acceptability by end users. Readability then becomes an important criteria. It is shown that a simple k-nearest-neighbor algorithm can achieve good performance with the help of the domain expert. As an illustration, a risk score made of only four attributes is presented for the french population.

I. INTRODUCTION

According to the World Health Organization, starting from 2010, cancer has become the leading cause of death worldwide [1] while, in western countries, new cases of cancer have grown at high rate in the last twenty years. Therefore, alongside clinical research to cure each cancer, prevention become a major concern in order to decrease its impact in our society. Among four existing levels of prevention, the general population is mostly interested in the two first levels. Primary prevention, that aims at avoiding occurrence of a disease and secondary prevention, that aims to diagnose and treat a disease as soon as possible, could both benefit from a widely used risk score.

Indeed, building a risk score for a given type of cancer can be a way to promote information and dialog between the medical profession and each patient [START_REF] Testard | The war on cancer[END_REF]. With appropriate information, anyone can improve its way of life by acting on specific risk factors to decrease its global risk. Early diagnosis is another way to sharply decrease death rate. Women with high risk could benefit from a quantified risk assessment that encourage them to follow screening program recommendations.

For example, women over 50 years old in France and over 40 years old in USA are recommended to perform a mammography every two years to detect early stage breast cancer which is the first cause of cancer for women. Generalized use of risk scores could benefit every patient, first by giving them solutions to lower their risk and second by convincing them to enter national prevention programs for specific diseases. Risk scores have to provide very good detection capabilities to identify high risk profiles among the general population. Risk scores also have to be built on readable modelling methods that can be easily understood by any patient or medical practitioner. Meeting our objectives with those constraints make the design of such a risk score challenging because it requires to deal with heavily imbalanced data (less than 2 % of breast cancer cases), using environmental factors that haven't high prediction power, with a restricted set of algorithms that meet readability criterion.

For breast cancer, the need for an efficient and readable risk score fit in a medical process currently under development by oncologists, epidemiologists and cancer treatment specialists. In this medical process, women should be able to check their own risk whether during a counseling appointment with a physician or a gynecologist, or through a dedicated internet website. Women with a high risk profile should then be encouraged to spend a day in a new generation risk clinic. There, a team of specialists will perform screenings for early stage breast cancer and determine her individual risk of developing a breast cancer by doing all tests and examinations in a single day to lower psychological impact. If needed, medical team may create a personalized plan to treat a detected cancer or to lower these risks.

This medical process is based on the ability to detect high risk french women that might need in depth counseling, but it lacks a readable and efficient risk score adjusted to french population. We suggest a method to conceive such risk score and a practical solution formalized as an information system adapted to data accessibility, computation capacity and deployment constraints. This article is organized in five sections. Section II provides a description of constraints of the medical process and software solutions. Section III describes existing risk models for breast cancer and source data we use. Section IV presents our data mining process and performances on the database. Conclusion is in section V.

II. FROM MEDICAL PROCESS TO INFORMATION SYSTEM

A. Medical process to improve prevention and screening First proofs of generalized early screening efficacity have been gained with breast cancer. Prevention treatments for high risk women have been proven to be efficient [START_REF] Fisher | Tamoxifen for prevention of breast cancer: Report of the national surgical adjuvant breast and bowel project p-1 study[END_REF]. Breast cancer research is also a field where treatments were highly improved with, as a result, a lower mortality rate. For this reasons, further screenings and prevention improvements could help to reduce mortality strongly. Goal of the medical process is to individualize both screening and prevention with help of a risk clinic, a risk clinic where high risk women will be advised to get an appointment for in depth analyze of their risks.

In order to identify women with higher risk in the general population, a risk score can be used. Two ways are considered to get women to have their risk assessed. First an internet website that any woman will be able to visit. Second, a software component integrated in the physician or gynecologist information system that provides alerts if a woman has an high risk profile depending on data available in her medical records. Both systems will advise user of the possibility to meet a team of specialists at a risk clinic.

Risk clinic have to be a decentralized structure that receive women in a care network of a given territory. Depending on the assessed risk for a woman, she may be directed to an appropriate care facility, a mammography center for example. For women with very high risks, a single appointment in a major care facility have to allow to: do a blood test for collection of clinical data, aggregate familial medical history and undergo basics screenings such as mammography. Depending on test results, the woman will get an individualized surveillance program or be sent to an appropriate care facility if a cancer is detected.

This medical process involves 6 stakeholders. The ordering party is the government that need more efficient processes to reduce mortality. The contractor is a non-profit institution in the medical field that will gather skills to ensure the ordering party that the several steps of the process will be coherent and functional. Epidemiologists (as a domain experts) and data miners have to design an efficient, easy to use and readable risk score. End users are the physician or gynecologist that will use the risk score, but also the woman that wants to assess her risk thanks to an internet website. Oncologists will lead in depth analyze for very high risk women whereas specialists will support prevention and screening for high risk profiles. At the end of the process, the beneficiaries will be the women who will get personalized advices or adapted treatments.

B. Constraints and impacts

The medical process described in II-A is based on the assessment of a women breast cancer risk. Practitioners keep total control over the advice given to the woman, but assistance can be provided through a risk score embedded in an information system that meets two constraint types.

Constraints on risk score: Risk score has to be efficient and readable to be used in a medical environment. Indeed, despite their performances, statistical breast cancer risk scores are not widely used in a prevention context. One of the reason may be that those scores are not understood by users. Users training is an explanation but all potential users will not be able to be trained in statistics. The risk score has to be based on an understandable and effective model. Readability and effectiveness impact the kind of algorithm we choose to build our risk score.

Risk scores have to allow inclusion of specific attributes for acceptability by end users. Physicians and specialists may have a priori ideas about good attributes of a model in a risk computation context. In the process of creating the risk score, an attribute selection step has to allow an expert, who has knowledge of users a priori ideas, to intervene on the choice of attributes.

Acquisition cost has to be very low. Breast cancer risk could be assessed with a better accuracy thanks to genetic risk factors (presence of BCRA1/2 gene for example) or with costly medical examination. But in front of a website or during a medical appointment, questions have to be both short and simple. Therefore, the risk score has to be based on environmental, reproductive or familial factors that are easy to answer.

Constraints on information system: List of attributes has to be flexible. Depending on the profile of the woman, some attributes may not be available or relevant to assess the risk. For example, the type of hormone replacement therapy may not be known at time of assessment. A risk score still have to be computed without the missing information. Another example is the age at menopause: if a women is not yet in menopause, age at menopause is not a relevant attribute to use, but still, risk score have to be computed. The system have to be able to deliver a risk score depending on the attributes that are available from user.

Delivering the score has to made almost in real time. For a statistically reliable risk score, model have to be based on a large enough learning set that increases computation time. But, to be able to deliver a risk level quickly, it has to be instantly available. As a consequence, the system has to precompute scores so that it can deliver a value without dealing with all the learning set each time a score is needed. Real time is a global constraint on the system, but it allows the learning step to be time consuming if needed, to increase efficiency or readability.

In epidemiology, data that allow to study cause of diseases are evolving permanently and are not publicly available. As a consequence, being able to add examples to the learning set as data are collected must be allowed by the system architecture, as well as not embedding a non-public database in a publicly available software. Having an online system that is interrogated by a client software through a web service is a solution to avoid distribution of the database and to allow availability of data.

To summarize, our information system has to provide a risk score that is efficient and readable, usable with missing values, available in real time while allowing to update the private database used to build it.

C. Information system overview

To meet constraints described in II-B that result from the medical process described in II-A, we propose different components built to be integrated in the system. They are designed to work with the information systems of public health departments and the information systems used by physicians in their offices, see Fig. 1.

In order to design and test risk scores, we need to ease access to:

• one or multiples database servers that stores different kinds of data about people who entered an epidemiological study. They contain scanned versions of paper questionnaires that were sent to people by mail, digitized versions of medical documents that were obtained from physicians and hospitals, database that contains raw answers obtained by optical character recognition on scanned questionnaires and digitalized medical documents. Database server also contains cleaned versions of a part of the data and generated data that were deducted from people answers. Access to these database servers is restricted because they contained sensitive data about people in a hospital environment. Access to these servers depends on the data format: raw data can be reached through a document managing system or a patient data management system and database tables can be accessed directly from a statistical software. • one or multiples computational servers. They have to provide enough power to compute all necessary data that make up our risk score and allow to run test procedures of the risk score. Useful data and executable programs have to be remotely sent on the server. Computation power can be reached through executable programs with a secured remote access. The data miner computer acts as a platform to access all information systems that provide needed resources under various forms. From this platform, data are gathered, explored, discretized and converted into a generic format. Modelling and risk score testing decisions are automated and implemented in a software. Such software is adapted to the computational server to take advantage of its multi-core and memory capabilities. Once a risk score model is chosen, it is made available through a database table that contains precomputed risk score levels. Indeed, to meet constraints of section II-B, we choose to allow its access through an internet solution.

Users (women, physicians, specialists) will be able to access the risk score through a web service from a web server. Access will be made easy: women will be able to calculate their risk score with a website whereas physicians or gynecologists will likely use it through their information systems thanks to components compatible with their medical software. To protect data, only precomputed risk score are uploaded on the web server that will provide a web service. It will be used to power a website for women and components that will be integrated into the different kinds of information systems owned by physicians.

III. BREAST CANCER: AVAILABLE RISK SCORES AND DATA

A. Breast cancer risk scores 1) Epidemiological approach: Thanks to the Breast Cancer Detection Demonstration Program, using an unconditional logistic regression model, Gail et al [START_REF] Gail | Projecting individualized probabilities of developing breast cancer for white females who are being examined annually[END_REF] conceived what is now the most commonly used model to predict breast cancer risk with environmental factors. Risk factors information was collected through a home interview from 6,000 women. Eventually, among 15 risk factor, only 5 were chosen to be part of the final model: age of patient, age at menarche (first natural menstrual period), age at first live birth, number of previous breast biopsies and number of first degree relatives affected by breast cancer. Cumulative risk of breast cancer is calculated by a multiplication of each of the five relative risks. Then, individual risk of a women can be computed by multiplying the cumulative risk by an adjusted population risk of breast cancer. The model was validated and adapted on several other population: asian and pacific islander american women [START_REF] Matsuno | Projecting individualized absolute invasive breast cancer risk in asian and pacific islander american women[END_REF], italian women [START_REF] Decarli | Gail model for prediction of absolute risk of invasive breast cancer: Independent evaluation in the florence-european prospective investigation into cancer and nutrition cohort[END_REF] and on american women [START_REF] Costantino | Validation studies for models projecting the risk of invasive and total breast cancer incidence[END_REF] of the CASH study (see below). On the 82,109 women of the Nurses Health Study, using Gail's model modified by Costantino [START_REF] Costantino | Validation studies for models projecting the risk of invasive and total breast cancer incidence[END_REF], Rockhill reports discrimination with an AUC (Area Under ROC Curve, see performance measurement in section IV-C) of 0.58 [START_REF] Rockhill | Validation of the gail et al. model of breast cancer risk prediction and implications for chemoprevention[END_REF].

Another well known risk score is the Claus model [START_REF] Claus | Autosomal dominant inheritance of early-onset breast cancer. implications for risk prediction[END_REF], also known as CASH model, based on the Cancer and Steroid Hormone Study lead by the Centers for Disease Control and Prevention. This genetic-based model was built using a segregation analysis for women who have a familial history of breast cancer. It aimed at understanding the breast cancer transmission model using familial data (mother and sister) of 9,418 women. Risk factors were: age of patient, number of first and second degree relatives affected by breast cancer. Individual risk can be easily obtained by reading a dedicated table depending on the number of affected relatives.

Using environmental factors, Barlow et al [START_REF] Barlow | Prospective breast cancer risk prediction model for women undergoing screening mammography[END_REF] proposed a risk prediction model based on the Breast Cancer Surveillance Consortium (BCSC) database which contains 2.4 millions screenings mammograms and answers to associated selfadministered questionnaires. Using two logistic regression models, a risk score was built with 4 or 10 risk factors depending on the menopausal status. Unlike the Claus model, it can be used for women without familial history of breast cancer. Compared to Gail's, it gains use of two attributes with higher prediction power: hormonal therapy and breast density. Barlow et al report AUC: 0.631 for premenopausal women based model and 0.624 for postmenopausal women.

2) Data mining approach: Because most approaches of breast cancer risk prediction deal with cancer relapse in the data mining field, authors did not faced as imbalanced data as epidemiological studies did (see III-A1). However, it is worth highlighting two significant studies involving mining algorithms and medical data.

Jerez-Aragonèz et al [START_REF] Jerez-Aragonés | A combined neural network and decision trees model for prognosis of breast cancer relapse[END_REF] chose the prognosis of breast cancer relapse to build a decision support tool. A database that gather information about 1,035 patients of the Malaga Hospital, Spain was used. Similar attributes than Gail (age, age at menarche and first full term pregnancy) were analyzed along biological tumor descriptors. Selection of most relevant prognosis factor was done with a tree induction based method. Attributes were then used to predict cancer relapse with an artificial neural network by computing a Bayes a posteriori probability in order to generate the prognosis system.

The Surveillance, Epidemiology and End Results (SEER) database was used by Endo et al [START_REF] Endo | Comparison of seven algorithms to predict breast cancer survival[END_REF] to implement common machine learning algorithms to predict survival rate for women affected by breast cancer. Even if positive examples were highly represented (18,5 %) in the database and area under the ROC curve wasn't chosen as performance metric, it is relevant to note that logistic regression had the highest accuracy, artificial neural network showed the highest specificity and J48 decision trees model had the best sensitivity.

Because we want to predict breast cancer risk for women among the general population, we face highly imbalanced data: new breast cancer cases rate is lower than 1 %. Usually, dealing with imbalanced data can be done both at algorithmic and data level according to Japkowicz et al [START_REF] Japkowicz | The class imbalance problem: A systematic study[END_REF] and Visa et al [START_REF] Visa | Issues in mining imbalanced data sets -a review paper[END_REF]. Indeed, guiding the data mining process or sampling data are solutions to increase detection performances for high risk profiles, but we will not use them in our approach due to the nature of the algorithm we chose.

Jerez-Aragonèz's [START_REF] Jerez-Aragonés | A combined neural network and decision trees model for prognosis of breast cancer relapse[END_REF] and Endo's [START_REF] Endo | Comparison of seven algorithms to predict breast cancer survival[END_REF] studies show how mining techniques can be used to build classification tools on medical databases while considering missing data and processes. But, even for epidemiological risk scores, they do not consider all of our constraints, such as readability, easy of use for patients, physicians and specialized doctors in their day to day interactions and adaptation to the french population.

B. Available data 1) E3N cohort:

To build a risk score for french women, we use data of one of the largest epidemiological cohort (a group of people with a shared characteristic) study in France: E3N (" Étude Épidémiologique des femmes de l' Éducation Nationale") [START_REF] Clavel-Chapelon | E3n, a french cohort study on cancer risk factors[END_REF]. It is a prospective study that includes women of the MGEN, "Mutuelle Générale de l' Éducation Nationale", a health insurance plan primarily covering employees of the french national education system. The E3N study is the french component of EPIC, the European Prospective Investigation into Cancer and Nutrition. Since 1990, 98,995 volunteers women, born from 1925 to 1950, were asked to fill 10 selfadministered questionnaires about their lifestyle (for example diet, reproductive factors, alcohol and tobacco consumption, physical activity, etc), regular use of medical treatments (hormonal treatments for example) and personal medical history As E3N is a prospective study, information was collected before any major pathology occurred, so women did not have cancer at time of inclusion in the cohort. When cancer cases are reported, invasive breast cancer cases are ascertained by obtaining histopathology reports. Up to 2011, invasive breast cancer cases were ascertained at a 92.4 % rate. Deaths information and causes of death, that may include unreported cancer cases, are obtained thanks to family members and the MGEN insurance database in accordance with the french data protection legislation.

2) Dataset construction: To conceive a breast cancer risk score that meets constraints of section II-B, we need a dataset with environmental data as attributes and cancer status as class. Several types of breast cancer exist, we define a women as being affected by breast cancer only if it is ascertained as an invasive breast cancer with no distinction depending on its hormone status. Among all attributes available from the E3N study, we decide to include in our dataset, only known risk factors. Risk factors are attributes for which we know from the literature, that epidemiologically speaking, they have an impact on the breast cancer risk. Impact should be understand at different levels: direct impact (e.g. hormone treatment), intermediate impact (e.g. age at menarche for hormone exposure), risk marker (e.g. number of biopsies). These 12 attributes are sum up in the Table I.

We decide to build a 5-years risk score in order to aggregate enough cancer cases and to offer a medium-term perspective to the woman and the physician. The dataset is built with attributes that describe a woman at time of answer of the fifth questionnaire in 1997: it is our baseline. As each women had several months to answer the 1997's questionnaire, each 5years timeframe does not start at the same time for each one of them. For each woman, cancer status is defined as positive if breast cancer occurs in a 5-years delay after her 1997 response, negative if not.

Dataset includes 92,078 women. During the 5-years timeframe, 1,647 women (1.79 %) have developed invasive breast cancer: they are labeled as positive examples whereas other 98.21 % are labeled as negative examples. Dataset is therefore strongly imbalanced. For some of the risk factors, questions were asked in previous questionnaires. Those questionnaires are used to retrieve information that is not available from the 1997's questionnaire: for example age at menarche is usually between 8 and 18 years old and it does not matter if information was collected 5 years before baseline, in 1992 instead of 1997. Remaining missing values, that could not be found in any questionnaire, are replaced by median value for the attribute. Discretization was realized in order to maintain distribution of values.

IV. BUILDING OUR RISK SCORE

A. CRISP-DM based process 1) Main objectives: The main objective of our approach is to build a method to conceive a risk score and a practical solution formalized as an information system that will be used through a business oriented software by a physician or through an internet website by any woman. As statistical models spread with difficulty in the medical field, we aim to find a model with good scoring performance and good readability. We say a model has a good readability if it can be understood thanks to a simple picture or sentence: it has to be quickly readable during a medical appointment in a primary prevention context. Furthermore, we have other constraints: patients need actionable attributes to change their lifestyle in a prevention context, physicians have a priori ideas about good attributes of a model in a risk computation context, both of them want immediately usable score (see II-B). In addition, a generic algorithm that can be easily adapted to various pathologies in different countries is desirable.

2) General process: Our approach follows the CRoss Industry Standard Process for Data Mining (CRISP-DM) [START_REF] Chapman | Crisp-dm 1.0 stepby-step data mining guide[END_REF] data-mining methodology. Figure 2 shows the 6 steps of this process where gray ones identify our major contributions.

Business understanding: An expert with knowledge of the needs of users help us to prioritize our objectives (see section IV-A1) and to assess the situation. We decide to focus on a scoring task (no classification or prediction).

Data understanding: The E3N cohort (see section III-B) contains numerous known breast cancer personal risk factors. Data are used by epidemiologists to discover associations between factor and cancer. Data preparation: To deal with highly imbalanced data, we can apply rebalancing algorithms on this data but it is not the focus of the paper. As our modeling step will be based on the negative/positive ratio, we want to minimize modification of data balance in order to provide risk score closest as possible as reality. The only modification we apply is normalization thanks to a division by the standard deviation. To limit overfitting, it was decided to apply a random 75/25% split between training and validation set.

Modeling: If several data mining algorithms were considered, domain expert suggested to use a k-nearest-neighbor algorithm because it uses a concept of similarity which is easily understandable by end-users without explaining a complex formula. Moreover, such algorithm is able to deal with imbalanced data if there is enough positive examples among neighbors. We generate models and search for the best combination of attributes by performing an exhaustive search (see section IV-B) on a limited set of combinations. The reason is that the expert issued a recommendation of using a restricted number of factors to make the risk score easy to use. Obviously, for large combinations, computation time can increase sharply, but it is not a problem as models are generated offline. As an example, for data described in section III-B1 and for 1 to 6 attributes combinations, it takes approximately 15 hours to compute.

The k-nearest-neighbor algorithm also meets the genericity criterion: as a non parametrical model, it fully relies on data (no statistical parameters optimization) so it can be adapted for other pathologies and other dataset from different countries by replacing the dataset.

Evaluation: We evaluate generated models with Receiver Operating Characteristic validation (see section IV-C) using Area Under Curve (AUC) in order to sorts models by scoring performance. Then, our expert has to choose the most useful models leveraging on the AUC performance combined with its knowledge of users needs. ROC evaluation of every generated model is automatized in our software as well as descriptive statistics to characterize the neighborhood: mean, median and standard deviation for both k value and number of positive examples. We are still improving our process to formalize and support expert choice through a dedicated graphical user interface. It provides useful information about each model and ease the browsing among hundreds of models by offering sorting and filtering options.

Deployment: We are currently working to develop web services that will power a web application and components for physicians information systems.

B. Focus on k-nearest-neighbor implementation

To provide experts with interesting models, k-nearestneighbor algorithm (see [START_REF] Fix | Discriminatory analysis, non-parametric discrimination: consistency properties[END_REF], [START_REF] Cover | Nearest neighbor pattern classification[END_REF]) is used with various size of attributes combinations (from 1 to 5 attributes) and several k values were used (see section V). Performance of each of hundreds generated combinations is tested for each values of k.

We implement the k-nearest-neighbor algorithm in two steps:

• Selection of neighborhood: for a combination of attributes (e.g. age and number of first degree relatives with breast cancer), a score value has to be computed for each combination of values (e.g. age=54 and number of relatives with breast cancer = 1). To compute such score value, a neighborhood has to be defined for each values combination. To determine if a profile of the database belong to the neighborhood of a combination of values, an euclidean distance is used to compute the distance between a combination of value and every single record of the dataset using a normalized version of the values. Thus, at least k of the nearest records of the database are included in the neighborhood. The neighborhood may not have always the same size because for a given group at the same distance, if k is not reached yet, all neighbors at the same distance are added to the neighborhood. • Scoring function: the score of a combination of values, is the ratio between the number of breast cancer cases (i.e. positive examples) and the size of the neighborhood. In epidemiology, the ratio of individuals having a disease in a population is called prevalence. This ratio was chosen because it is well known by physicians, easily explainable to a risk score user and it is directly built on the number of patient diagnosed with breast cancer among patients with a similar profile.

C. Focus on ROC evaluation

Our performance metric has to depict how positive instances are assigned with higher scores than negative ones: we used the Receiver Operating Characteristic (ROC) [START_REF] Egan | Signal detection theory and ROC analysis[END_REF] to measure performance due to the continuous nature of our classifier. The ROC curve enable to measure detection performances using a moving threshold to classify examples of the validation set. Moreover, it allows direct comparison with epidemiologicalbased scores from the literature.

ROC mechanism: negative examples labeled as positive by the algorithm are called a false positives whereas positive examples labeled as positives are called true positives. The ROC curve is plotted with the false positive rate on the X axis and the true positive rate on the Y axis [START_REF] Fawcett | An introduction to ROC analysis[END_REF], both rates are calculated for a given threshold. ROC curve can be summarized in one number: the Area Under the ROC Curve (AUC). The area is a portion of the unit square, its value is in an [0,1] interval. The best classifier will have an AUC of 1.0 (i.e. all positive examples are assigned with higher score than negative ones) whereas an AUC of 0.5 is equivalent to random score assignment.

Each k value of each tested combination of attributes is assigned with a ROC curve and the corresponding AUC in order to help the expert to choose the best model.

V. EXPERIMENTAL RESULTS

On a publicly available database [START_REF] Barlow | Prospective breast cancer risk prediction model for women undergoing screening mammography[END_REF], we have empirically shown that a data mining approach can provide better performances from a discrimination and a readability perspective in [START_REF] Gauthier | Breast cancer risk score: a data mining approach to improve readability[END_REF]. In this section, we report the results obtained with a dataset built on the E3N cohort with specific objective: build a risk score for french women that will be used in a global medical process called a risk clinic.

A. Expert knowledge limits dataset

Before testing performances of our k-nearest-neighbor (knn) implementation, we submit our list of 12 attributes (see section III-B2) to a domain expert. With his knowledge of physicians and specialists a priori ideas about a good composition for a breast cancer risk score (a constraint of section II-B), he advices us to keep only 8 attributes among 12 (see Table II). He chooses this 8 attributes, or a subset of them, because they are commonly used by physicians to roughly assess their patients risk and because their are recognized in the community to be good risk factor for breast cancer assessment. The first list of attributes was objectively conceived with knowledge of impact on cancer risk as only criterion. It did not depend on constraints set by process stakeholders. However, this second list is created to meet a specific constraint of the process. It is build subjectively to take into consideration, the knowledge that we have from users for the specific purpose of building a breast cancer risk for a defined set of users. We call this 8 attributes dataset, the restricted dataset.

B. Scoring performances

An experiment set was designed to test how the knn algorithm performs on the restricted dataset build from the E3N cohort. As one of our constraints is to build a readable risk score, we decide to limit the size of the combination to 5 attributes. We select all combinations with a size s of 1 to 5 attributes among n = 8 attributes of the restricted dataset, meaning we have Among combinations of one attribute (see Table IV), the hrt combination is the best factor to predict breast cancer with an AUC of 0.572. Next best attribute is age with an AUC of 0.552. Usually, age is the best predictor for breast cancer risk but, as the hormone replacement treatment is given after the menopause, hrt attributes carry a partial information about age of the woman plus the kind of hormone replacement treatment used by women of the cohort. Epidemiologically speaking, domain expert states that this combination of one attribute has no sense and performance is not good enough, so we consider other available combinations.

Among tested combinations, the best performance is achieved by the hrt, age, kdeg1, fftp, agemeno combination. The attribute only slightly increase performances, then to meet our readability constraint, we choose not to use it. Same reasoning could apply to the 4-attributes combination compared to the 3-attributes combination, but the 3-attributes combination does not display use of age. Acceptability constraint would not be met without age attribute according to our

C. Details on the chosen combination

Combination hrt, age, kdeg1, fftp is chosen because it meets our readability and acceptability constraints with a high level of performances compared to other combinations. Figure 3 shows ROC curve for the combination.

Neighborhood characteristics: To assign a risk score to examples of the validation set, all scores were computed on the learning set first by defining a neighborhood for each profile (the profile age=52, hrt=0, fftp=23 and kdeg1=1 for example), and then by computing the prevalence value (see section IV-B) for the profile. For this chosen combination, the generated neighborhoods contain at least 9,000 neighbors and a maximum 12,359 of neighbors (because k is not a strict value, see section IV-B). On average, it contains 9,296.7 neighbors (median: 9,161, standard deviation: 373,2). As the score is based on the number of breast cancer cases among neighbors, it is interesting to look at the number of cancer cases in the neighborhoods: at least 110 and a maximum of 288 cancer cases in the neighborhood, on average 235.0 cases (median: 251 and standard deviation: 38.0).

Performance stability: In order to run a knn algorithm, the size of neighborhood has to be set. Since only k closest neighbors are used to compute the ratio healthy vs. diseased, risk score value depends on k value. If the neighborhood is too small, few breast cancer cases are included and if the neighborhood is too large, patient profiles are too different: in both cases the risk score is not reliable. Figure 4 presents evolution of the AUC depending on the size of the neighborhood.

With an undersized neighborhood, performances are low but then, as k increases, performances increase with a maximum of 0.604. From 6,900 to 10,300 neighbors, performances are always higher than 0.600 meaning that the algorithm p with p = 3 to 6 instead of an euclidian distance with p = 2, has not provided better AUC measures on the chosen combination. Adding the distance-weighted knn rule described by Dudani [START_REF] Dudani | The distance-weighted k-nearest-neighbor rule[END_REF] do not increase performances with the weighting function that were tested on the chosen combination. We decide not to try further as we do not want to overlearn.

D. Discussion

Our intention is to build an IS that integrate: access to cohort database and a computational server to conceive of a risk score that will be used to detect women with high risk for breast cancer. To achieve this purpose, a major step is to build a risk score that fit in the following medical process: risk assessment, cancer related advices, risk clinic appointment.

By testing every combination of attributes from a restricted dataset build with help of a domain expert and by carefully selecting a set of attributes depending on its discrimination capacity and the nature of the risk factors to meet our constraints, we choose the hrt, age, kdeg1, fftp combination. With an AUC of 0.604, the risk score have similar performances than other scores that use similar risk factors as Rockhill et al [START_REF] Rockhill | Validation of the gail et al. model of breast cancer risk prediction and implications for chemoprevention[END_REF] with Gail's model on an american population. But performances are not as high as those obtained by Barlow et al [START_REF] Barlow | Prospective breast cancer risk prediction model for women undergoing screening mammography[END_REF] or by Gauthier et al [START_REF] Gauthier | Breast cancer risk score: a data mining approach to improve readability[END_REF] on the american Breast Cancer Surveillance Consortium database. Domain expert explains this difference by the absence of attributes with higher prediction power such as breast density. This information being currently retrieved from women of the E3N cohort, we hope to increase performances of our risk score.

Discrimination performances may be increased using another knn enhancements than those we tested, but we might prefer to focus on improving another aspects of our IS. For example, use of the software to generate combinations and test them could benefit from tracability and reproductibility improvements in order to keep track of dataset, knn settings and associated results. Graphical user interface for end users also need work to ease the use of the risk score and make the knn concept as readable as possible.

Nevertheless, our study has some limitations. First, even if we used one of the few databases large enough to be representative of the targeted population, findings from a database based on employees of the national education system require cautious extrapolation to general population. Conversely, as our risk score is build on a ratio between number of cancer cases and size of neighborhood, it depicts association between a woman profile ant a breast cancer risk. Over-or underrepresentation of profile in the database, compared to general population, should have limited impact on extrapolation of the built risk score. Second, use of expert knowledge could be improved: process followed by domain expert to select useful combinations could benefit from adapted graphical user interface to browse output of our knn implementation.

VI. CONCLUSION

In order to create an efficient and readable breast cancer risk score for french women that fit in a medical process built to detect high risk women that may benefit from in depth counseling and regular screenings, we studied some constraints that have to be considered to design an information system to build risk scores for major pathologies. We then instantiate the system with a knn algorithm and the domain expert knowledge.. Discrimination performances of the knn algorithm for the combination chosen by domain expert are in the same range than measures reported by studies on different population with similar attributes thanks to widely used statistical methods. Moreover our proposal meets deployment and readability constraints that make the challenge to produce a practical and readable solution usable by physicians and women.

Fig. 1 .

 1 Fig. 1. Applicative architecture of information system

  n-s)! = 218 combinations to test. Each one of the 218 combinations was tested with 39 values of k nearest neighbors, meaning 8,502 configurations were computed.

Fig. 3 .

 3 Fig. 3. ROC curve for the hrt, age, kdeg1, fftp combination

  TableIIIprovides descriptive statistics by combination size. Each attributes combination appears only once in the statistics with the best AUC computed among all values of k that were tested. As the size of combinations increases, AUC increases more slowly. Decline of the standard deviation shows that combination tends to have more similar performances because combinations tends to include more risk factors with prediction power.

			TABLE III
		DESCRIPTIVE STATISTICS BY COMBINATION SIZE
	Size Combinations AUC Mean AUC Std Dev. AUC Median
	1	8	0.533	0.023	0.530
	2	28	0.556	0.022	0.553
	3	56	0.569	0.019	0.565
	4	70	0.577	0.017	0.579
	5	56	0.584	0.016	0.589

ACKNOWLEDGMENT

The authors are indebted to all participants for providing the data used in the E3N study and to practitioners for providing pathology reports. They are grateful to Mrs. R. Chaït, Mrs. M. Fangon, Mrs. L. Hoang, Mrs. C. Kernaleguen and Mrs. M. Niravong for their technical assistance.


Cédric Boutillier
email: boutilli@ccr.jussieu.fr

Non colliding paths in the honeycomb dimer model and the Dyson process

In this paper we describe a natural family of random non-intersecting discrete paths in the dimer model on the honeycomb lattice. We show that when the dimer model is going to freeze, this family of paths, after a proper rescaling, converge to the extended sine process, obtained traditionally as the limit of the Dyson model when the number of particles goes to infinity.

Introduction

The Dyson model is a dynamic version of the Gaussian Unitary Ensemble in random matrix theory. It models a collection of identical charged particles diffusing, but repelling one another by electrostatic interaction. The trajectories of these particles are thus forced not to intersect. If we put these particles in an harmonic well of potential, the system reaches an equilibrium. When the number of particles grows to infinity, and in equilibrium, the system is invariant by translation in time and in space. It is described by a determinantal process, known as the extended sine process.

The dimer model belongs to discrete statistical mechanics. This system was introduced in the 1930's [START_REF] Fowler | Statistical theory of perfect solutions[END_REF] to model the adsorption of diatomic molecules on the surface of a crystal, and consists in a probabilistic study of perfect matchings of a graph. It is remarkable that this model is exactly solvable when no supplementary interaction between edges is imposed. The partition function as well as correlation functions can be expressed explicitly, and have a determinantal form: this is a model of discrete free fermions. We will look at a particular case, when the graph is (a subgraph of) the honeycomb lattice G.

It has been noticed several times that dimer models have some properties very similar to those of random matrices. In presence of boundary conditions, the dimer model may exhibit several behaviours. In some regions (called liquid regions, usually far from the boundary) correlations decay polynomially, and in some cases, the kernel used to 1 compute correlations between dimers along a given line is the discrete sine kernel (see for example [START_REF] Baik | Uniform asymptotics for polynomials orthogonal with respect to a general class of discrete weights and universality results for associated ensembles: announcement of results[END_REF][START_REF] Johansson | Non-intersecting paths, random tilings and random matrices[END_REF][START_REF] Okounkov | Random skew plane partitions and the pearcey process[END_REF]). Under a proper scaling, the continuous determinantal sine process is recovered, describing the bulk of large Gaussian Hermitian matrices [START_REF] Johansson | Non-intersecting paths, random tilings and random matrices[END_REF][START_REF] Boutillier | The bead model & limit behaviors of dimer models[END_REF]. In other regions, the configuration is frozen: there the configuration is almost deterministic. The behaviour of the system near the frontier between a frozen and a liquid region are related to the Airy process and thus to the statistics of the largest eigenvalue of a large Gaussian Hermitian matrix [START_REF]The arctic circle boundary and the Airy process[END_REF][START_REF] Ferrari | Step fluctuations for a faceted crystal[END_REF][START_REF] Okounkov | Random skew plane partitions and the pearcey process[END_REF]. See also [START_REF] Johansson | Eigenvalues of gue minors[END_REF][START_REF]The birth of a random matrix[END_REF] for other interesting results enforcing the relations between these objects. That is why dimer models are sometimes considered as discrete analogues of random matrices [START_REF] Kenyon | What is . . . a dimer?[END_REF].

In this article, we add to that list another example by relating the dimer model on the infinite honeycomb lattice with the Dyson model. A dimer configuration on the honeycomb lattice G can be described by the paths obtained by following the non horizontal dimers. The dimer model on G is endowed with a two-parameter family of ergodic Gibbs measure that are conditionally uniform. The two parameters can be seen as the drift of a path, and the average vertical distance between two successive paths.

We show that the extended sine process, obtained at the limit of the Dyson process in the bulk and in the steady regime, when the number of particles becomes infinite, is also obtained as the limit of Gibbs measures on dimer configurations of the infinite honeycomb lattice, after a proper scaling, when the drift is zero and the average distance goes to infinity.

Theorem 1. The family of discrete non-colliding paths in the honeycomb dimer model converges weakly, after a proper scaling, to the extended sine process.

The precise meaning of this convergence is explained in the beginning of Section 3. A possible strategy is to use Lindström-Gessel-Viennot [START_REF] Gessel | Binomial determinants, paths, and hook length formulae[END_REF]'s approach (see also [START_REF] Karlin | Coincidence probabilities[END_REF]) to describe the combinatorics of the non-intersecting paths in terms of determinants. However, we will make use of the technology introduced by Kasteleyn [START_REF] Kasteleyn | Dimer statistics and phase transitions[END_REF] and developed in [START_REF] Kenyon | Local statistics of lattice dimers[END_REF][START_REF] Kenyon | Dimers and amoebae[END_REF] to compute correlations between dimers.

The Dyson model, and thus the extended sine process can be obtained as the scaling limit of another statistical mechanical model, called the random vicious walks introduced by Fisher in [START_REF] Fisher | Walks, walls, wetting, and melting[END_REF] (see for example [START_REF] Katori | Functional central limit theorems for vicious walkers[END_REF]). In that model, a finite number of walkers on Z start at time T 0 from even integers and end at T 1 . The measure on these paths is the uniform measure on all configurations without intersections. What differs in this work from the vicious walkers case is that due to the ergodicity property of the Gibbs measures on dimer configurations, the family of path is already at the discrete level infinite, and invariant by (discrete) translations.

The paper is structured as follows: in Section 1, we give a brief description of the Dyson model, and the extended sine process, its limit when the number of particles goes to infinity, in the steady regime. In Section 2, we present some aspects of the dimer model on the infinite honeycomb lattice and its Gibbs measures, as well as the correspondence with a family of discrete non-colliding paths. In Section 3, we prove the weak convergence of the probability measures on the random paths in the dimer model to the extended sine process. The proof of two propositions used in that section are postponed to Section 4 and Section 5.

The Dyson model and the sine process

The Dyson model [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF] is defined as the motion of n identically charged diffusive particles λ 1 (t), . . . , λ n (t) repelled one from another by an electrostatic force, but attracted to the origin by an elastic force. This is a Markov process whose evolution kernel p(t, λ, µ) satisfies the following equation

∂p ∂t = 1 2 k j=1 ∂ 2 p ∂λ 2 j - ∂ ∂λ j ∂Φ(λ) ∂λ j p , (1)
where the potential

Φ(λ) = -n j=1 λ 2 j + i<j log |λ j -λ i |. When n = 1, this is the Ornstein-Uhlenbeck process.
This process can be interpreted as the evolution of the eigenvalues of a random Hermitian matrix, where the diagonal entries evolve as an Ornstein-Uhlenbeck process with speed 2 and the real and imaginary parts of the entries above the diagonal , as an Ornstein-Uhlenbeck process with speed 1, all these motion being independent. The statistics of the eigenvalues at a fixed time are those of a random matrix from the Gaussian Unitary Ensemble (GUE) [START_REF] Mehta | Random matrices[END_REF]. This process is an example of a determinantal process where the correlation functions are given by determinants of the extended Hermite kernel, a space-time extension of the Hermite kernel appearing in the analysis of the GUE.

If we scale the Dyson process in the bulk when n goes to infinity, we get the sine process, a determinantal process associated with the extended sine kernel defined as

S(t -t, y -y) =        1 0 e -(t -t)φ 2 /2 cos φ(y -y) dφ π if t ≤ t - +∞ 1 e -(t -t)φ 2 /2 cos φ(y -y) dφ π if t > t (2)
When t = t , this reduces to the usual sine kernel sin(y -y) π(y -y) . The sine process describes an infinite system of particles repulsing each other by Coulomb force, and which is homogeneous both in space and time. Such a system has been studied in [START_REF] Dyson | A Brownian-motion model for the eigenvalues of a random matrix[END_REF][START_REF] Spohn | Interacting Brownian particles: a study of Dyson's model[END_REF][START_REF] Katori | Infinite systems of non-colliding Brownian particles[END_REF].

Statistics of occupation of the process can be expressed in terms of S: given τ 1 < • • • < τ k and (I j ) k j=1 a sequence of Borel sets of R, the probability that the number N I j (τ j ) of curves crossing I j equals n j at time τ j , for all j ∈ {1, . . . , k} is given by the following expression

P ∀ j ∈ {1, . . . , k}, N I j (τ j ) = n j = p∈N k (-1) |p| p! n! • • • (τ 1 ,I 1 ) n 1 +p 1 ו••×(τ k ,I k ) n k +p k det |n+p| S(ξ i -ξ j ) d n+p ξ (3)
where ξ i stands for the variables (τ, y) ranging over the ith factor of the Cartesian product (τ 1 , I 1 )

n 1 +p 1 × • • • × (τ k , I k ) n k +p k .
In this expression, we used multi-index notations:

n! = k j=1 n j !, |n| = k j=1 n j .
It will be convenient in the sequel to use also for z = (z 1 , . . . , z k )

z n = z n 1 1 • • • z n k k , ∂ n ∂z n = ∂ n 1 ∂z n 1 1 • • • ∂ n k ∂z n k k ,
and

{τ } × I = {τ 1 } × I 1 , . . . , {τ k } × I k ), ({τ } × I) n = ({τ 1 } × I 1 ) n 1 × • • • × ({τ k } × I k ) n k .
The probabilities in (3) can be seen as derivatives of a generating function for the joint distribution of N I 1 (τ 1 ), . . . , N I k (τ k ) : posing

Q(z) = Q {τ }×I (z) = E   k j=1 (1 -z j ) N I j (τ j )   = n∈N k (-z) n n! E N I j (τ j )! N I j (τ j ) -n j )! = n∈N k (-z) n n! ({τ }×I) n det |n| S(ξ i -ξ j ) d n ξ,
we have

P ∀ j ∈ {1, . . . , k}, N I j (τ j ) = n j = (-1) |n| n! ∂ n ∂z n Q {τ }×I (z) z=(1,...,1) . ( 4
)
Thus, all the information about the N I j (τ j ) is encoded in the function Q {τ }×I .

The dimer model on the honeycomb lattice

Let G be the honeycomb lattice. Since all the cycles are even, there is a bipartite coloring of the vertices of G in white and black. The group Z 2 acts by color-preserving translations on G, generated by x and ŷ, and the fundamental domain contains one white vertex w 0 and one black vertex b 0 . They are represented in Figure 1. All the vertices of G are images of vertices of the fundamental domain under some translation of Z 2 . We label them by their color and element of Z 2 . We thus denote by

w x,y = w 0 + (x, y), b x,y = b 0 + (x, y).
the white vertex and the black vertex in the fundamental domain with coordinates (x, y). The edges of G are classified in 3 types, a, b and c, according to their orientations around a white vertex (see Figure 1b).

A dimer configuration C of G is a subset of edges such that every vertex of G is incident with exactly one edge of C. An example of a dimer configuration on a piece of the honeycomb lattice is presented in Figure 2. We now define probability measures on the set of all dimer configurations of G.

Gibbs measures and correlations in the honeycomb dimer model

When the graph is finite, one can easily define the uniform probability measure on dimer configurations. But since the honeycomb lattice is an infinite periodic graph, it is not possible to construct directly a uniform measure. However, it was proved in [START_REF] Kenyon | Dimers and amoebae[END_REF][START_REF] Sheffield | Random Surfaces: Large Deviations Principles and Gradient Gibbs Measure Classifications[END_REF] that there exists a two-parameter family of Gibbs measures on the set of dimer configurations, having the following properties:

• they are ergodic under the action of Z 2 by translation,

• they are conditionally uniform: if the configuration is fixed in some annular region of the graph, the configurations outside and inside the annulus are independent, and the measure induced inside is the uniform measure.

These measures are obtained as weak limits of grand-canonical probability measures on dimer configurations on tori G/nZ 2 , with different activities a, b and c for the different kinds of edges. The measure P a,b,c obtained this way depends only on the ratios b/a and c/a. This model is exactly solvable. The free energy as well as the correlation functions can be computed explicitly. Following Kasteleyn's ideas [START_REF] Kasteleyn | Dimer statistics and phase transitions[END_REF][START_REF]Graph theory and crystal physics[END_REF], Kenyon [START_REF] Kenyon | Local statistics of lattice dimers[END_REF] proved that under P a,b,c , the probability that the edges e 1 = (w

x 1 ,y 1 , b x 1 ,y 1 ), . . . , e n = (w x n ,y n , b x n ,y n )
are covered by dimers in the random configuration has a determinantal form and is given by the following formula

P a,b,c e 1 , • • • , e n = n i=1 weight(e i ) det 1≤i,j≤n K -1 (b x i ,y i , w x j ,y j ) . ( 5
)
where weight(e) equals a, b or c depending of the type of the edge e, and

K -1 (b x ,y , w x,y ) = K -1 (b x -x,y -y , w 0 ) = |z|=|w|=1 z -(y -y) w x -x a + 1 w (b + cz) dz 2iπz dw 2iπw . (6)
K -1 is the inverse of the weighted bipartite adjacency matrix K, whose rows are indexed by white vertices and columns by black vertices.

K w,b = weight(e) if e = (w, b), 0 otherwise.
Kasteleyn discovered that this matrix K encodes all the combinatorics of the model. K is named after him the Kasteleyn operator. K -1 is called the inverse Kasteleyn operator.

The entries of K -1 are computed using Fourier transform where z (resp. w) is the multiplicator associated with translations in the y direction (resp. -x).

The decay of the entries of K -1 is directly related to the presence of singularities of the integrand on the unit torus T 2 = {(z, w) ∈ C 2 ; |z| = |w| = 1}. Depending on the weights, the behavior of the integral defining K -1 , and thus properties of the model can be drastically different:

• If a weight is greater than the sum of the two others, then with probability 1, we only see edges of one type (the type with highest weight). We say that the system is frozen. There are three frozen phases, each one corresponding to a type of edges.

• If (a, b, c) satisfy the triangle inequalities, then the density of each type of dimers is positive with probability 1, and correlations between dimers decay polynomially with the distance. The system is said to be in the liquid phase.

We will be interested in the phase transition between the liquid phase and the frozen a-phase, when b = c and a approaches b + c from below.

Non colliding paths in the honeycomb dimer model

The dual graph of the honeycomb lattice is the triangular lattice. This duality allows to define a bijection between dimer configurations of the honeycomb lattice and tilings of the plane by rhombi with a π 3 angle: one can associate to each edge the rhombus obtained by the union of the two dual faces corresponding to the endpoints of the edge. It will be convenient to use those two points of view.

The three types of rhombi, corresponding to the three types of edges, are also called a, b and c. To describe a tiling, it is sufficient to give the positions of tiles of type, for Figure 2: An example of a dimer configuration on a piece of the honeycomb lattice G. The corresponding tiling with rhombi is drawn in dotted lines.

example, b and c, forming left-to-right non-intersecting paths1 (see Fig. 3). A measure on dimer configurations can be viewed equivalently as a measure on tilings on the plane as well as a measure on non-colliding left-to-right paths.

When a approaches b + c from below, one sees more and more a-tiles, and the b/cpaths are more and more dilute : the average distance between two paths goes to infinity. The probability of seeing a piece of such a path in a finite window goes to zero as we approach the phase transition. To see something non trivial in the limit, one needs to rescale correctly the lattice in the same time as we let a go to b + c. To find the correct regime, it can be enlightening to look at the behavior of one single path.

Looking at a single path near the transition

In what follows, we will take b = c = 1 and a = 2 cos ε 2 . We will denote by P ε the Gibbs measure for these weights and K -1 ε the corresponding inverse Kasteleyn operator. For this particular choice of weights, after the change of variables u = w/ √ z and e iθ = z, K -1 ε has the following form:

K -1 ε (b x,y , w 0 ) = K -1 ε (x, y) = = π -π e iθ(-y+x/2) S 1 u x 2 cos(ε/2)u + 2 cos(θ/2) du 2iπ dθ 2π
.

The pole at u = -cos(θ/2) cos(ε/2) is outside of the unit disk if and only if 2 cos(θ/2) > a, i.e. if and only if θ ∈ (-ε, ε). The integral over u in the definition of K -1 ε can be computed Figure 3: A piece of a tiling of the plane with rhombi, in bijection with a dimer configuration on G. One can distinguish two left-to-right paths made of rhombi of type b and c corresponding to non horizontal dimers.

using residues, and we get the following expressions for the entries of that operator

K -1 ε (x, y) =        (-1) x+1 2 cos(ε/2) [-ε,ε] cos(θ/2) cos(ε/2) x e iθ(-y+x/2) dθ 2π if x < 0 (-1) x 2 cos(ε/2) [-π,π]\[-ε,ε] cos(θ/2) cos(ε/2) x e iθ(-y+x/2) dθ 2π if x ≥ 0 (7)
From ( 5), the probability of seeing an edge of type a is aK -1 ε (0, 0) = 1 -ε π . As a goes to 2, the non-intersecting paths become more and more distant, and one can believe that in the limit, these paths will behave almost like simple random walks, at least locally. This is indeed the case, as stated in the following lemma: Lemma 1. When b = c = 1, and in the limit a → 2 -, the left-to-right paths locally behave like simple random walks. More precisely, for any fixed n > 0, the distribution of n successive steps -up ( ) or down ( ) -in a left-to-right path converges to the uniform measure on { , } n . Proof. Let (0, y 0 ),. . . ,(n, y n ) be a sequence of Z 2 such that |y k+1 -y k | = 1 for all k ∈ {0, . . . , n -1}. For the sake of simplicity, we will write b j and w j instead of b j,y j and w j,b j . The sequence (y k+1 -y k ) can be interpreted as a sequence ω = (ω n ) of n steps: +1 for up, and -1 for down.

The fact that a left-to-right path of rhombi covers b 0 and then follows the sequence of steps ω can be formulated in terms of dimers as follows: the edge (w 0 , b 0 ) is not covered by a dimer, and the edges (b 0 , w 1 ),. . . ,(b n-1 , w n ) are covered by dimers.

Thus by formula (5), the probability that, given that a path covers b 0 , it follows the succession of steps ω, can be expressed as a determinant:

1 1 -aK -1 ε (b 0 , w 0 ) det    K -1 ε (b 0 , w 1 ) • • • K -1 ε (b n-1 , w 1 ) . . . . . . . . . K -1 ε (b 0 , w n ) • • • K -1 ε (b n-1 , w n )    . ( 8
)
Using the explicit expressions for K -1 ε given in [START_REF] Fowler | Statistical theory of perfect solutions[END_REF], one sees that

K -1 ε (b j , w j ) = 1 2 + o(1), and K -1 ε (b j , w i ) = ε 2 + o(ε) for j < i.
The determinant in ( 8) is equal up to higher order terms to det

        0 1 2 • • • 0 0 1 2 . . . . . . . . . . . . . . . 0 0 1 2 ε 2 0 • • • • • • 0         = ε 2 n
and hence, the probability (8) equals

1 ε + o(ε) ε 2 n + o(ε) = 1 2 n + o(1) (9)
When ε goes to zero, this conditional probability converges to the corresponding probability for a simple random walk.

Close to the phase transition, we are left with an infinite collection of random paths that individually behave locally like a random walk, but that are conditioned not to intersect with one another.

This individual random walk behavior is an indication that if we want to get a nontrivial model in the limit, we have to choose the same regime as the one giving a non trivial limit for the simple random walk: the diffusive regime.

Convergence to the extended sine process

For a fixed ε, we represent our paths of rhombi by continuous functions as follows. We first embed the honeycomb lattice H in the plane correctly. Recall that the vectors defining the coordinates on H are not orthogonal and that the horizontal direction is that of the vector ŷ -2x. In order to rescale the horizontal direction by a factor ε 2 and the vertical direction by a factor ε, we map every white vertex with coordinates (x, y) to the point with usual cartesian coordinates ε 2 8 x, ε 2 (2y -x) . This diffusive scaling seems natural since the paths of rhombi behave like random walks at small scales.

A path of rhombi is then encoded by the piecewise linear path joining all the white vertices covered by this path. We get a bi-infinite family of paths (X ε n (t)) n∈Z , indexed by time t ∈ R. By convention, the path corresponding to the index 0 is the one passing the closest to y = 0 at time 0.

Such a family of paths is in bijection with a random tiling with rhombi. The Gibbs measure P ε on tilings (or equivalently on dimer configurations of G) can therefore be viewed as a probability measure on C(R ∞ ), the set of continuous functions on R with values in

R ∞ = x = (x j ) j∈Z |∀j ∈ Z x j ∈ R . The space R ∞ is a complete separable metric space with the metric dist(x, ȳ) = +∞ k=-∞ 1 2 |k| |x k -y k | 1 + |x k -y k |
The topology induced by this distance is the topology of coordinatewise convergence, and the sets

B k,ε (x) = ȳ ; ∀ j ∈ {-k, . . . , k} |y j -x j | < ε
form a basis for this topology. One can then define a topology on C(R ∞ ), using the || • || ∞ norm with respect to that distance on R ∞ , which allows to define the notion of weak convergence of probability measures on C(R ∞ ): a sequence of probability measure (P ε ) on C(R ∞ ) converges weakly to P if for all functions f on C(R ∞ ) continuous for the topology defined, lim

ε→0 E ε [f (X)] = E[f (X)],
where E ε and E are the expectations with respect to P ε and P respectively, and X represents the random family of paths. The aim of this section is to prove the main result of this paper Theorem 2. The probability measures (P ε ), as measures on bi-infinite families of continuous paths on R, converge weakly in C(R ∞ ) to the sine process.

The proof of the convergence of the ensemble of paths to the sine model is based on precise asymptotics of the inverse Kasteleyn operator when ε goes to zero, that are proved in Lemma 2. These asymptotics allow us to prove that the finite dimensional distributions of P ε converge to those for the sine process in Proposition 1. Then we prove that (P ε ) is tight in Proposition 2. A key point of the proof of this proposition is a simple comparison between probability for paths of rhombi and for simple random walks given by Lemma 3.

Since the proofs of Proposition 1 and 2 are rather long and technical, they are exposed separately, in Section 4 and 5 respectively.

Let us begin with the asymptotics of the kernel in the scaling regime described above, when ε goes to zero.

Lemma 2. For every

ε > 0, let (x ε , y ε ) ∈ Z 2 such that lim ε→0 1 8 ε 2 x ε = T ∈ R, lim ε→0 1 2 ε(2y ε -x ε ) = Y ∈ R. ( 10
)
Recall that a = 2 cos ε 2 .

• If T = 0, then lim ε→0 (-1) xε+1 ε aK -1 ε (x ε , y ε ) = e T S(T, Y ). ( 11
) • If x ε = 0 and y ε = 0, lim ε→0 -a ε K -1 ε (0, y ε ) = sin(Y ) πY = e 0 S(0, Y ). ( 12
)
• If x ε = 0 and y ε = 0,

lim ε→0 1 -aK -1 ε (0, 0) ε = 1 π = e 0 S(0, 0). ( 13
)
Moreover, the quantities aK -

1 ε (x ε ,y ε ) ε (or 1-aK -1 ε (0,0) ε if (x ε , y ε ) = (0, 0)) are uniformly bounded in ε.
This lemma shows that the discrete kernel for the dimer model on the honeycomb lattice converges, up to a factor e T and except for the third case, to the extended sine kernel. It is based on the exact formulae for K -1 ε given in Equation [START_REF] Fowler | Statistical theory of perfect solutions[END_REF].

Proof. Suppose first that T < 0. In this case, for ε small enough, x ε < 0 and

K -1 ε (x ε , y ε ) = (-1) x ε+1 a ε -ε 2 cos(θ/2) a xε e i(-y ε +x ε /2)θ dθ 2π = (-1) x ε +1 a ε -ε cos(θ/2) cos(ε/2) 8T ε 2 (1+o(1)) e -iθ Y ε (1+o(1)) dθ 2π = ε(-1) x ε +1 a 1 -1 cos εφ 2 cos ε 2 8T ε 2 (1+o(1)) e -iY φ(1+o(1)) dφ π = ε(-1) x ε +1 a 1 -1 1 - ε 2 (φ 2 -1) 8 + o(ε 2 ) 8T ε 2 (1+o(1))
e -iY φ+o (1) dφ 2π

When ε goes to zero, the integrand converges to e T e -φ 2 T e -iY φ . A simple application of Lebesgue dominated convergence theorem shows then that

lim ε→0 (-1) x+1 ε aK -1 ε (x ε , y ε ) = e T 1 -1 e -φ 2 T e -iY φ dφ 2π = e T 1 0 e -φ 2 T cos(Y φ) dφ π = S(T, Y ). ( 14
)
When T > 0, then x ε > 0 for ε small enough, and we have

K -1 ε (x ε , y ε ) = (-1) x ε a [-π,π]\[-ε,ε] 2 cos(θ/2) a x ε e i(-yε+xε/2)θ dθ 2π (15) = (-1) x ε a π ε 2 cos(θ/2) a xε cos (-y + x/2)θ dθ π . ( 16
) ( 17
)
If we change the variable in the integral, defining φ = εθ, we get

K -1 ε (x ε , y ε ) = ε(-1) x ε a π/ε 1 cos εφ 2 cos ε 2 8T ε 2 (1+o(1)) cos(Y φ + o(1)) dφ π = ε(-1) x ε a A 1 cos εφ 2 cos ε 2 8T ε 2 (1+o(1)) cos(Y φ + o(1)) dφ π + ε(-1) x ε a π/ε A cos εφ 2 cos ε 2 8T ε 2 (1+o(1)) cos(Y φ + o(1)) dφ π
for a given large number A. The second integral can be bounded by

ε(-1) xε a π/ε A cos εφ 2 cos ε 2 8T ε 2 (1+o(1)) cos(Y φ + o(1)) dφ π ≤ Ce -T (A-1)+o(1) (18)
uniformly in ε. Lebesgue's dominated convergence theorem shows that the first integral converges to

e T A 1 e -φ 2 T cos(Y φ) dφ π ( 19
)
As a consequence, letting A go to ∞, we get

lim ε→0 (-1) x ε +1 ε aK -1 ε (y, x) = -e T +∞ 1 e -φ 2 T cos(Y φ) dφ π . ( 20
)
When x ε = 0, K -1 (x ε , y ε ) given by

K -1 ε (0, y ε ) = 1 a π ε cos(y ε θ) dθ π . ( 21
)
If

y ε = 0, then -a ε K -1 ε (0, y ε ) = - a ε sin(πy ε ) -sin(εy ε ) aπ = sin(Y ) πY + o(1). ( 22
)
However, if

y ε = 0, K -1 ε (0, 0) = π-ε aπ and thus lim ε→0 1 -aK -1 ε (0, 0) ε = 1 π . ( 23
)
The uniform bounds are obtained directly form the expression of K -1 ε given by ( 7). When T < 0,

a ε K -1 ε (x ε , y ε ) = 1 ε ε -ε 2 cos(θ/2) a x ε e iθ(-y ε +x ε /2) dθ 2π ≤ 1 π 1 -ε 2 2 -x ε /2 ≤ C. ( 24
)
When T > 0,

a ε K -1 ε (x ε , y ε ) = 1 ε [-π,π]\[-ε,ε] 2 cos(θ/2) a x ε e iθ(-y ε +x ε /2) dθ 2π ≤ 1 π 1 -ε 2 2 -x ε /2 π ε 1 e -ε 2 x ε 8 t 2 dt ≤ C. ( 25
)
When x ε = 0, the bound is obtained by a direct evaluation of [START_REF] Fowler | Statistical theory of perfect solutions[END_REF].

As in the case of the extended sine process, for given times τ 1 , . . . , τ k and unions of intervals I 1 , . . . , I k , we can compute for every ε the probability that for every j, the number N ε I j (τ j ) of paths in I j at time τ j equals n j . The following proposition states that these probabilities converges to the corresponding probabilities for the sine process.

Proposition 1. The finite dimensional distributions of X ε converge to that of the extended sine process.

For all τ 1 , . . . , τ m ∈ R, and I 1 , . . . , I m ⊂ R finite unions of intervals,

lim ε→0 P ε N ε I 1 (τ 1 ) = n 1 , . . . , N ε I k (τ k ) = n k = P [N I 1 (τ 1 ) = n 1 , . . . , N I k (τ k ) = n k ] (26)
We prove also tightness of the family of probability measures (P ε ) in C(R ∞ ). Since R ∞ is a nice complete separable metric space, tightness is characterized by the following two properties [START_REF] Billingsley | Convergence of probability measures[END_REF] we have to check:

Proposition 2.
i) The sequence of distributions of (X ε (0)) is tight: for each η > 0, there exists a sequence of closed intervals

([-b k , b k ]) k∈Z such that ∀ε, P ε [∃ k ∈ Z, X ε k (0) / ∈ [-b k , b k ]] < η (27)
ii) For each L > 0, for each positive δ and η, there exists an α ∈ (0, L) such that

∀ε, P ε    sup -L/2<s,t<L/2 |s-t|<α dist X ε (s), X ε (t) > δ    < η (28)
Tightness and convergence of finite dimensional distributions of the family (P ε ) imply the weak convergence of this family to a probability measure on C(R ∞ ) whose finite distributions coincide with that of the sine process. The family of discrete random paths therefore converges weakly to the family of random continuous paths described by the sine process. This ends the proof of Theorem 2, modulo the proof of the two last propositions.

Proof of proposition 1

This section is devoted to the proof of the convergence of finite dimensional distributions of the family (P ε ). This is done by studying the convergence of generating functions Q ε for the distribution of the number of visits of a finite number of intervals at given times when ε goes to zero.

We denote by τ ε the quantity τ 8ε 2 and (I) ε τ the set of integers

(I) ε τ = y ∈ Z |ε(y - 1 2 τ ε ) ∈ I .
The set of white vertices with coordinates (x, y) ∈ {τ ε } × (I) ε τ are the closest to {τ } × I.

The function

Q ε of the variable z = (z 1 , . . . , z k ) ∈ C k defined by Q ε (z) = E ε   k j=1 1 -z j ) N ε I j (τ j )   = n∈N k E ε   k j=1 N ε I j (τ j )! N ε I j (τ j ) -n j !   (-z) n n! ( 29
)
where

z n = k j=1 z n j
j , is a generating function for the probabilities we are interested in. Indeed, from (4), we have:

P ε N ε I 1 (τ 1 ) = n 1 , . . . , N ε I k (τ k ) = n k = p∈N k (-1) p n! p! E ε   k j=1 N ε I j (τ j )! N ε I j (τ j ) -(n j + p j ) !   (30) = (-1) n n! ∂ n ∂z n Q ε (z) z=(1,...,1) . ( 31
)
We will show that Q ε converges to the corresponding generating function for the extended sine process. For this, we first compute the quantities E ε k j=1

N ε I j (τ j )! N ε I j (τ j )-n j !
and their limit when ε goes to 0. These quantities are given by the sum over all distinct n j -tuples of white vertices in {τ ε j } × (I j ) ε τ j for j ∈ {1, . . . , k} of the probability that these white vertices are covered by a path of rhombi of type b and c, or equivalently, that the horizontal edges with weight a incident with these vertices are not present in the random dimer configuration. We have

E ε   k j=1 N ε I j (τ j )! N ε I j (τ j ) -n j !   = j=1,...,k y (j,1) ,...,y (j,n j ) ∈(I j ) ε τ j distinct P ε a τ ε 1 ,y (1,1) , . . . , a τ ε 1 ,y (1,n 1 ) , . . . , a τ ε k ,y (k,n k ) / ∈ C
(32) where a x,y is the horizontal edge with weight a incident with the white vertex with coordinates (x, y) in H, and C is the random dimer configuration.

Using multilinearity of the determinant and an inclusion-exclusion argument, one can write each of these probabilities

P ε a τ ε 1 ,y (1,1) , . . . , a τ ε 1 ,y (1,n 1 ) , . . . , a τ ε k ,y (k,n k ) / ∈ C (33) as a determinant det δ i,j -aK -1 ε (τ ε i -τ ε j , y i -y j ) = det |n|          1-aK -1 ε (τ ε 1 -τ ε 1 ,y (1,1) -y (1,1) ) ••• -aK -1 ε (τ ε k -τ ε 1 ,y (k,n k ) -y (1,1) ) . . . . . . . . . . . . . . . . . . . . . -aK -1 ε (τ ε 1 -τ ε k ,y (1,n 1 ) -y (k,n k ) ) ••• 1-aK -1 ε (τ ε k -τ ε k ,y (k,n k ) -y (k,n k ) )          (34)
where i and j belong to the following list of indices (1, 1), . . . , (1, n 1 ), . . . , (k, 1), . . . , (k, n k ) , (

and the angular brackets i of a couple i = (i 1 , i 2 ) represent its first factor i 1 .

The determinant is unchanged if we multiply each column i by e -τ i (-1)

τ ε i and each line by e τ j (-1) -τ ε j . It follows from Lemma 2 that this determinant is asymptotic to ε |n| det S(τ i -τ j , Y ε i -Y ε j ) , (36)
where

Y ε j = ε(y j -τ ε j ) ∈ I j .
The sum (32) becomes a Riemann sum converging as ε goes to zero to the following integral

• • • (τ 1 ,I 1 ) n 1 ו••×(τ k ,I k ) n k det S(ξ i -ξ j ))d n ξ. ( 37
)
Therefore, the coefficients of Q ε (z) converge to those of Q(z), the generating function for the sine process. By Lemma 2, the entries of the matrix in (34) are uniformly bounded, say by M . By applying Hadamard's lemma , we bound the sum (32) by

E ε   k j=1 N ε I j (τ j )! (N ε I j (τ j ) -n j )!   ≤   k j=1 I j n j   (|n|M 2 ) |n|/2 (38)
uniformly in ε. Thanks to the additional factor 1 n! , the general term of the series in (30) is uniformly bounded for ε sufficently small and z in a compact set, by the general term of an absolutely convergent series. Thus, by Lebesgue's dominated convergence theorem, all the derivatives of the series Q ε (z) converge uniformly on compact sets of C k to the corresponding derivatives of the generating function Q(z) for the sine process. In particular,

lim ε→0 P ε N ε I 1 (τ 1 ) = n 1 , . . . , N ε I k (τ k ) = n k = P [N I 1 (τ 1 ) = n 1 , . . . , N I k (τ k ) = n k ]
where P is the probability measure for the extended sine process, with the kernel S given by (2).

Proof of proposition 2

This section is devoted to the proof of the tightness of the family of probability measures (P ε ) ε>0 . First, we check the point i) of proposition 2. Let η > 0. We look for a sequence of segments such that

∀ε, P ε [∃ k ∈ Z, X ε k (0) / ∈ [-b k , b k ]] < η
is satisfied. If we take a sequence verifying that

∀ k ∈ N b k = b -k and b k+1 > b k (39)
then the condition ( 27) is equivalent to

∀ε > 0 P ε [∃ k ∈ N, N b k < 2k + 1] < η (40)
where

N b k = N ε [-b k ,b k ] (0)
is the number of paths crossing the interval [-b k , b k ] at time 0. We will use Chebyshev inequality to get a bound on the probability (40). For this, we need first to estimate and bound the average and the variance of the random variable

N L = N ε [-L,L] (0
) for any L > 0. The average of N L under P ε is given by

E ε [N L ] = y∈[-L ε , L ε ]∩Z P ε [there is a path crossing at (0, εy)] (41) = y∈[-L ε , L ε ]∩Z (1 -P ε [a y ∈ C]), ( 42
)
where a y = a 0,y is the horizontal edge whose white end has coordinates (0, y) in H. There are 2 L ε + 1 integer points in [-L ε , L ε ], and

P ε [a y ∈ C] = 1 -ε π . Therefore, E ε [N L ] = 2 L ε + 1 ε π = 2L π 1 + o(1) . ( 43
)
Moreover, the variance of N L is given by

Var ε (N L ) =E ε N 2 L -E ε [N L ] 2 (44) =E ε [N L (N L -1)] -E ε [N L ] 2 + E ε [N L ] < E ε [N L ] (45) Indeed, E ε [N L (N L -1)] -E ε [N L ] 2 ,
given by the following expression

y,y ∈[-L ε , L ε ]∩Z y =y P ε a y , a y / ∈ C -    y∈[-L ε , L ε ]∩Z P ε [a y / ∈ C]    2 = (46) = y,y ∈[-L ε , L ε ]∩Z y =y det 1-aK -1 ε (0,0) -aK -1 ε (0,y-y ) -aK -1 ε (0,y -y) 1-aK -1 ε (0,0) - y∈[-L ε , L ε ]∩Z 1-aK -1 ε (0,0) 2 (47) = - y,y ∈[-L ε , L ε ]∩Z y =y a 2 K -1 (0, y -y ) 2 -(1 -aK -1 ε (0, 0)) 2 2 L ε + 1 < 0 (48)
is a sum of negative terms. By Chebyshev inequality, we get

P ε N L < L π ≤ P ε N L -E ε [N L ] < L π ≤ Var ε (N L ) (L/π) 2 ≤ 2π L 1 + o(1) . (49) Thus, for η ≤ 3π 2 2 , taking b k = 2π 3 (k+1) 2 3η
, using the fact that for k ≥ 0,

2k + 1 ≤ (k + 1) 2 ≤ b k π ,
we finally get that

P ε [∃ k ∈ N, N b k < 2k + 1] ≤ +∞ k=0 P ε [N b k < 2k + 1] (50) ≤ +∞ k=0 P ε N b k < b k π (51) ≤ +∞ k=0 2π b k (1 + o(1)) = η 2 (1 + o(1)) ≤ η, (52)
for ε small enough, uniformly in η. Therefore, the family of distributions of (X(0)) under P ε is tight. The proof of point ii) goes as follows. First, as the topology on R ∞ is that of coordinatewise convergence, the condition (28) is equivalent to

∀ L > 0, ∀ η > 0, ∀ k ∈ N, ∃ α ∈ (0, L), ∀ ε > 0, ∀ j ∈ {-k, . . . , k} P ε     sup s,t∈(-L 2 , L 2 ) |s-t|<α |X j (s) -X j (t)| > δ     < η. (53)
It is therefore sufficient to check this condition for each path separately. It is well-known that this property is true for simple random walks: for any δ > 0, there exists an α > 0 such that the number of properly rescaled walks (Y ε ) with N = L/ε steps such that sup

s,t∈(-L 2 , L 2 ) |s-t|<α |Y (s) -Y (t)| > δ (54)
is less than η2 N uniformly in ε. Therefore, by the following, rather technical lemma 3, since the probability of a sequence of steps of X j is bounded by a constant times the same probability for a simple random walk, we get by decomposing the event realisation by realisation of the path X j from 0 to L,

P ε     sup s,t∈(-L 2 , L 2 ) |s-t|<α |X j (s) -X j (t)| > δ     < C L • η, ( 55
)
where the constant C L depends only in L.

Lemma 3. Let L > 0. There exists a c L > 0 such that for every ε small enough and every (ω n ) 0≤n≤ L/ε sequence of upward and downward steps, the probability that the first L/ε steps of a path of rhombi X ε j , conditioned on its starting position (x 0 , y 0 ), coincide with (ω n ), is bounded by

P ε first L ε steps of X ε j coincide with (ω n ) X ε j starts at (x 0 , y 0 ) ≤ c L 2 L/ε . (56)
If the number of steps was finite, one could have proceeded as in Lemma 1. However, here we want a control for a number of steps going to infinity as ε goes to zero.

Proof. Define g(x, y) = 1 2 ε -ε cos θ 2 x e iθ(-y+x/2) dθ 2π . ( 57
) Since π -π 2 cos(θ/2) x e iθ(-y+x/2) dθ 2π = x k=0 x k π -π e i θ 2 (2k-x) e i θ 2 (-2y+x) dθ 2π = x y , (58)
K -1 (x, y) can be rewritten as

K -1 ε (x, y) = -2 a x+1 g(x, y) if x < 0, -2 a x+1 g(x, y) -1 2 x+1 x y if x ≥ 0. ( 59
)
Let N = L ε . Suppose that the sequence (ω n ) corresponds to a path covering the vertices

b 0 = b (0) y 0 , w 1 = w (1) y 1 , b 1 = b (1) y 1 , . . . , b N -1 = b (N -1) y N -1 , w N = w (N ) y N , (60)
with y i+1 -y i = 1+ω i 2 . The probability that a path goes through these vertices is given by

P ε [ there exists a path through b 0 , w 1 , . . . , b N -1 , w N ] = det K -1 ε (b j-1 , w i ) = det K -1 ε ((j -1) -i, y j-1 -y i ) (61)
where the indices i and j refer to white and black vertices respectively and range from 1 to N . Multiplying simultaneously each column i of the matrix (61) by ( -a 2 ) -i and each row j by ( -a 2 ) j does not affect the determinant. This manipulation cancels the coefficients -2 a (j-1)-i+1 in the expression of K -1 (b j-1 , w i ), leading to a simpler expression for the determinant (61).

P ε [there is path through b 0 , w 1 , . . . , b N -1 , w N ] = det g((j -1) -i, y j-1 -y i ) - δ j>i 2 j-i j -1 -i y j-1 -y i , (62)
where δ j>i = 1 is j > i and 0 otherwise. The matrix, the determinant of which we want to compute, is the sum of a matrix with entries g(x, y) = O(ε), and a strictly upper triangular matrix.

We have now to find an upper bound for the value of this determinant to get the estimate stated in the lemma. We will use the Hadamard inequality, but a direct application of it would give too approximative a bound. We will have to make some more manipulations on the rows and columns of this matrix.

First, note that for every x ∈ {-N, . . . , N }, and y = y or y + 1,

uniformly in x. Therefore, replacing for j ∈ {2, . . . , N }, column C j by C j -C j-1 , we get on the diagonal and under it entries that are O(ε 3 ) (except in the first column). It is important to observe that after these operations, the module of the entries (i, j) strictly above the diagonal does not increase, at least at the leading order in ε, and are bounded from above by a constant times (j -i) -1/2 , since

1 2 x+1 x y - 1 2 x x -1 y ≤ 1 2 x+1 x y ≤ 1 2 x+1 x x/2 ≤ 1 √ 2πx . ( 64
)
Now, we use the second column to put 0 instead of the first entries of columns C j , j > 2 by making the substitution C j ← C j + α j,2 C 2 for a suitable value of α j,2 in the determinant. From the bound on coefficients, we see that |α j,2 | = O (j -2) -1/2 and that therefore, the module of the entries of column C j , j ≤ 3 does not increase more than by O (j -2) -1/2 ε 3 = o(1).

(65)

We then continue this procedure, and for 3 ≤ k ≤ N -1, we use column C k to eliminate the entries of the k -1th row on columns C j , j > k.

After this succession of operations, the entries (i, j) of the matrix with j > i + 1 are 0. The module of the other entries of column C j , j ≤ 2 increased at most by

O ε 3 √ j -2 + O ε 3 √ j -3 + • • • + O ε 3 j -(j -1) = O(ε 3 j) ( 66
)
The probability we are interested is given by the determinant det

         O(ε) 1 2 + O( √ 2ε 3 ) 0 • • • 0 O(ε) O( √ 2ε 3 )
√ N ε 3 ) O(ε) O( √ 2ε 3 ) • • • O( √ N ε 3 )          . ( 67
)
Hadamard inequality states that a determinant is bounded by the product of the 2 -norm of its rows. In this particular case, it gives (71) Since the probability that there is a path of rhombi starting from b 0 is of order ε, we obtain the wanted bound for the conditional probability P ε the first L ε steps of X ε j coincide with (ω n ) X ε j starts at (x 0 , y 0 ) ≤

c L 2 N . ( 72
)
for some constant c L depending on L. This ends the proof of the last lemma, giving the last argument to end the proof of theorem 2.

Figure 1 :

 1 Figure 1: a) A piece of the honeycomb lattice. The fundamental domain is represented by the shaded region. The two vectors x and ŷ generate the group of translations preserving the graph G and its bipartite coloring. b) Weights a, b, c are assigned to edges according to their orientation around a white vertex.

g(x + 1 ,θ 2 ) x e i θ 2 (-2y+x) e ±i θ 2 cos θ 2

 1222 y ) -g(x, y)

2 . 2 + O(ε 2 )

 222 P ε [there is path through b 0 , w 1 , . . . , b N -1 , w N ] (68) Recall that N = L ε = O 1 ε .Thus, from the fact that, for j ∈ {1, . . . , N -1}, finallyP ε [there is path through b 0 , w 1 , . . . , b N -1 , w N ] ≤ O(ε) × 1

These non intersecting paths can be seen directly in the dimer picture by superimposing on top of the random configuration, the frozen a-configuration.

Acknowledgments: We are grateful to Richard Kenyon for enlightening discussions on dimer models, which have motivated this work.