Juliette Stehlé

Nicolas Voirin

Alain Barrat
email: alain.barrat@cpt.univ-mrs.fr

Ciro Cattuto

Lorenzo Isella

Jean-François Pinton

Marco Quaggiotto

Wouter Van Den Broeck

Corinne Régis

Bruno Lina

Philippe Vanhems

High-resolution measurements of face-to-face contact patterns in a primary school Authors

Background

Little quantitative information is available on the mixing patterns of children in school environments. Describing and understanding contacts between children at school would help quantify the transmission opportunities of respiratory infections and identify situations within schools where the risk of transmission is higher. We report on measurements carried out in a French school (6-12 years children), where we collected data on the time-resolved face-to-face proximity of children and teachers using a proximity-sensing infrastructure based on radio frequency identification devices.

Methods and Findings

Data on face-to-face interactions were collected on Thursday, October 1 st and Friday, October 2 nd 2009. We recorded 77,602 contact events between 242 individuals (232 children and 10 teachers). In this setting, each child has on average 323 contacts per day with 47 other children, leading to an average daily interaction time of 176 minutes. Most contacts are brief, but long contacts are also observed. Contacts occur mostly within each class, and each child spends on average three times more time in contact with classmates than with children of other classes. We describe the temporal evolution of the contact network and the trajectories followed by the children in the school, which constrain the contact patterns. We determine an exposure matrix aimed at informing mathematical models. This matrix exhibits a class and age structure which is very different from the homogeneous mixing hypothesis.

Conclusions

We report on important properties of the contact patterns between school children that are relevant for modeling the propagation of diseases and for evaluating control measures. We discuss public health implications related to the management of schools in case of epidemics and pandemics. Our results can help define a prioritization of control measures based on preventive measures, case isolation, classes and school closures, that could reduce the disruption to education during epidemics.

Introduction

The role of children in the community spread of respiratory infections such as influenza is a challenging epidemiological issue [START_REF] Longini | Estimating household and community transmission parameters for influenza[END_REF][START_REF] Viboud | Risk factors of influenza transmission in households[END_REF]. Children are thought to act as a reservoir as they are the first to be infected by various respiratory infections such as influenza, respiratory syncytial virus, or rhinovirus. They can be infected at school because of the numerous close contacts occurring between school children, and then act as sources of infection into their households from where infections can further spread in the community [START_REF] Longini | Estimating household and community transmission parameters for influenza[END_REF][START_REF] Viboud | Risk factors of influenza transmission in households[END_REF]. However, little is known about how children actually mix in a school environment [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF][START_REF] Mikolajczyk | Social contacts of school children and the transmission of respiratory-spread pathogens[END_REF][START_REF] Glass | Social contact networks for the spread of pandemic influenza in children and teenagers[END_REF]. An accurate description and understanding of the contacts between children at the school level would help to quantify the transmission opportunities of respiratory infections, and to identify the situations during school days where the risk of transmission is higher. Ultimately, one goal of such quantification would be to assess which control and containment measures have the best performance. In addition the availability of quantitative descriptions of the contact patterns between young individuals has the potential to inform mathematical models that aim at describing the propagation of diseases in populations.

In order to reduce this knowledge gap, the research priorities comprise collecting data on activities and interactions of children, in particular in schools. Until recently, most empirical studies have relied on self-reported information such as questionnaire-based surveys to determine mixing patterns [START_REF] Zagheni | Using time-use data to parametrize models for the spread of close-contact infectious diseases[END_REF][START_REF] Wallinga | Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents[END_REF][START_REF] Mossong | Social contacts and mixing patterns relevant to the spread of infectious diseases[END_REF][START_REF] Beutels | Social mixing patterns for transmission models of close contact infections: exploring self-evaluation and diarybased data collection through a web-based interface[END_REF][START_REF] Hens | Mining social mixing patterns for infectious disease models based on a two-day population survey in Belgium[END_REF]. Such surveys may however suffer from biases due to self-reporting.

Recent advances in distributed sensing systems, based on mobile devices and wearable sensors, provide new ways of gathering data on human contacts and allow to mine the proximity relations and close-range interactions of individuals in real-world large-scale settings [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF][START_REF] Hui | Pocket switched networks and human mobility in conference environments[END_REF][START_REF] Eagle | Reality mining: sensing complex social systems[END_REF][START_REF] O'neill | Instrumenting the city: Developing methods for observing and understanding the digital cityscape[END_REF][START_REF] Pentland | Honest Signals: how they shape our world[END_REF][START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF][16][START_REF] Alani | Live Social Semantics, 8th International Semantic Web Conference ISWC[END_REF]. In this framework, a recent study [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF] has for instance given important insights into the proximity patterns of students, teachers and staff in a US high school. As the data available to a broad research community remain scarce, additional measurement campaigns in different settings and covering different schools, countries, and age groups are much needed, in particular to test for common patterns and differences, and to understand their public health implications.

We deployed a proximity-sensing infrastructure based on radio-frequency identification devices (RFID) in a French primary school, and used it to collect, in an unsupervised manner [START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF], time-resolved data on the face-to-face proximity of children and teachers. We report the number and duration of contacts, and their dependence on age, class, daytime and school spatial structure. Based on the results, we highlight specific situations where children are in contact and during which infections may be transmitted, and discuss the potential impact of containment measures such as class or school closures during seasonal epidemics or pandemics of respiratory infections.

Methods

Setting

The study took place in a primary school in Lyon, France during two days in October 2009. The age of the students (elementary cycle, composed of 5 grades) ranges between 6 and 12 years. In this school, each of the 5 grades is divided in two classes, for a total of 10 classes. Each class has an assigned room and an assigned teacher. The smallest class has 22 children and the largest 26, for a total of 241 children and 10 teachers. 232 children and all teachers participated in the data collection. The school day runs from 8.30am to 4.30pm, with a lunch break from 12pm to 2pm, and two breaks of 20-25 min around 10.30am and 3.30pm. Lunches are served in a common canteen, and a shared playground is located outside the main building. As the playground and the canteen do not have enough capacity to host all the students at the same time, only two or three classes have breaks at the same time, and lunches are taken in two consecutive turns.

Study design/Ethics statement and privacy

The French national bodies responsible for ethics and privacy, the Comission Nationale de l'Informatique et des Libertés (CNIL, http://www.cnil.fr) and the "Comité de Protection des personnes" (http://www.cppsudest2.com/), were notified of the study, which was approved by the relevant academic authorities (by the 'directeur de l'enseignement catholique du diocese de Lyon', as the school in which the study took place is a private catholic school). In preparation for the study, parents and teachers were invited to a meeting in which the details and the aims of the study were illustrated. Verbal informed consent was then obtained from parents, teachers and from the director of the school. All participants were given a Radio-Frequency Identification (RFID) badge and were asked to wear it at all times. Special care was paid to the privacy and data protection aspects of the study: The communication between RFID badges, the readers, and the computer system used to collect data were fully encrypted. No personal information of participants was associated with the identifier of the corresponding RFID badge. The only piece of information associated with the unique identifier of the badge was the class the corresponding individual was associated with.

Data collection infrastructure

The measurement infrastructure, developed in the context of the SocioPatterns project [16], is based on active RFID devices, embedded in unobtrusive wearable badges. Detailed information on how this technology is used to monitor social interactions and to identify contact patterns is available in Refs. [START_REF] Alani | Live Social Semantics, 8th International Semantic Web Conference ISWC[END_REF][START_REF] Isella | What's in a crowd? Analysis of face-to-face behavioral networks[END_REF][START_REF] Isella | Close encounters in a pediatric ward: measuring face-to-face proximity and mixing patterns with wearable sensors[END_REF]. Individuals are asked to wear the devices on their chests, so that badges can exchange radio packets only when the individuals wearing them face each other at close range (about 1 to 1.5 m). This range was chosen as a proxy of a close-range encounter during which a communicable disease infection can be transmitted, for example, either by cough or sneeze, or directly by hands contact. The infrastructure parameters are tuned so that the proximity of two individuals wearing the RFID badges can be assessed with a probability in excess of 99% over an interval of 20 seconds [START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF]16]. This time scale allows an adequate description of person-to-person interactions that includes brief encounters. We define that a "contact" occurs between two individuals during an interval of 20s if and only if the RFID devices worn by the individuals exchanged at least one radio packet during that interval. After a contact is established, it is considered ongoing as long as the devices continue to exchange at least one such packet for every subsequent 20s interval. Conversely, a contact is considered broken if a 20s interval elapses with no exchange of radio packets. The detected proximity relations are relayed from the RFID devices to receivers installed in the monitored environment and eventually stored in a central system..

Collected data

Data on face-to-face interactions between 232 children (96% coverage) and 10 teachers (100% coverage) across 10 classes were collected over two days (Thursday, October 1 st 2009 and Friday, October 2 nd 2009, see Table 1). Data were collected from 8.45am to 5.20pm on the first day, and from 8.30am to 5.05pm on the second day. Contacts were not recorded outside of these time intervals. Radio receivers (RFID readers) covered all the classrooms, the canteen, the stairways, and the playground. No information on contacts taking place outside the school or during sports activities was gathered.

Data analysis

The patterns of contacts between children, and the corresponding mixing patterns between classes and age groups are analyzed through several quantities describing the number of contacts between individuals, the duration of these contacts, and the cumulated time spent in contact by each pair of individuals, as well as their statistical distributions characterized in particular by average and coefficient of variation squared ( 2CV , defined as the squared ratio of the standard deviation to the mean of the distribution). A More precisely, we define the following weights quantifying the proximity relations of a pair of individuals i and j:

the occurrence ij o is equal to 1 if at least one contact event between i and j was recorded during the measurement period, and 0 otherwise; this is a symmetric quantity ( ).

As children are grouped into classes, we also compute:

the total number of contacts between children of classes A and B, [START_REF] Mossong | Social contacts and mixing patterns relevant to the spread of infectious diseases[END_REF]. All these quantities are computed both for the entire study duration and separately for each day.

∑ ∈ ∈ ij B j A, i AB n = n (for B A ≠ ; for B A = we have ∑ ∈A j i, ij AA n = n 2
In order to study the temporal structure of the interaction network, we also measure for each child i, and as a function of time, the number i k of other distinct children whom s/he was in contact with, as well as the total time i s the child i spent in contact with other children since the beginning of the measurements. In other terms, i k is the degree of node i in the contact network aggregated since the start of the first day, and i s is the corresponding node strength, i.e., the sum of the weights of the links inciding on i. In each case, we distinguish between contacts with children in the same class (

in i k , in i s
) and with children of different classes (

out i k , out i s ).
In terms of the proximity relations defined above, these quantities are therefore defined according to the following formulae (we consider in these formulae the example of a child i of class A):

∑ j ij i o = k , ∑ j ij i w = s , ∑ ∈A j ij in i o = k , ∑ ∈A j ij in i w = s , ∑ ∉A j ij out i o = k , ∑ ∉A j ij out i w = s
, where ij o and ij w are computed since the beginning of the first day. The temporal evolution of these quantities indicates how fast the overall contact patterns build up between the children, both at the class and at the school level.

We also build contact networks aggregated on a 20-minute timescale: each day is divided in sliding windows of 20 minutes, starting at intervals of 5 minutes and, for each 20-minute period, edges are drawn between those pairs of individuals for which at least one contact was recorded during this period. The average degree of each 20-minute network gives the average number of individuals with whom a given individual has been in contact with during the corresponding time window. By using these 20-minute sliding windows we filter out the fast fluctuations of the dynamical contact network and only retain the slowly-varying information on the network evolution.

As a summary of the contacts of each day, we additionally build two daily aggregated networks in which edges are drawn between a pair of individuals whenever at least one contact was recorded for that pair during the considered day. Each edge is weighted by the total time the corresponding individuals spent in contact during that day.

The two daily aggregated networks are compared using various measures. We compute the Pearson correlation coefficients between the characteristic parameters (number of contacts, total time spent in contact, etc.) measured for each individual in day 1 and in day 2. We also compare the network structures at a more detailed level, measuring the similarity between the neighborhoods of each node across the two days. A simple measure of similarity is given by the respective numbers of new and repeated distinct persons contacted in day 2 with respect to day 1.

This can be further refined by specifying if these new and repeated contacts occur within the same class or with individuals of other classes. Moreover, as each link i-j in the aggregated network is weighted by the total time i and j spent in face-to-face proximity (denoted by 1 , ij w for day 1 and by 2 , ij w for day 2), the similarity between the neighborhoods of an individual i in days 1 and 2 can be also quantified by the cosine similarity

( ) ∑ ∑ ∑ j ij j ij j ij ij w w w w = i sim 2 2 , 2 1 , 2 , 1 , ) ( ) (
This quantity is 1 if i had contacts in both days with exactly the same individuals, spending the same fraction of time in proximity of each one, and 0 if, on the contrary, i had contact in day 2 with totally different persons with respect to day 1.

Finally, by measuring the relative rates at which the RFID readers receive the packets emitted by individual badges, it is possible to perform approximate localization of the badges, and tell which RFID reader is closest to any given badge. Since the readers were installed in the classrooms, in the canteen, and in the courtyard, it is possible to detect in which of these areas each badge was situated at any point in time. This allows to construct the trajectories that children followed in space as they move within the school premises.

Results

Number of contacts

We recorded a total of 77,602 contact events involving 242 individuals (37,414 contacts on day 1 and 40,188 on day 2), with an average of about 317 contacts per individual on the first day (coefficient of variation squared

Duration of contacts

Most contacts are of short duration, but contacts of very different durations are observed, including rather long ones. Figure 3 shows the distribution (number of events in each bin divided by the bin width) of the contact durations P(t), defined as the fraction of contacts with duration t. Figure 3 highlights that large fluctuations are observed in the contact durations. While the average duration of a contact is 33 seconds, and 88% of the contacts last less than 1 minute, more than 0.2% of the contacts exceed 5 minutes. Overall, no characteristic contact time scale can be defined (the squared coefficient of variation of the distribution is

1 . 1 2 = CV ).
The heterogeneity of contact patterns is also observed for cumulated contact durations. shown by the whitish diagonal. Two-by-two light blocks around the diagonal also show that larger number and durations of contacts are observed between children of the same grade rather than with other grades. Finally, a separation between smaller grades (1 st to 3 rd ) and upper grades (4 th and 5 th ) grades is also apparent, most probably due to the lunch break schedule. The numerical values corresponding to Figures 4 and5 are given in Table 2, while Table 3 reports the average daily numbers and durations of contacts of an individual of one grade with individuals of another grade. Values are higher on the diagonal, the separation between smaller and upper grades is apparent, and values tend to decrease when moving away from the diagonal. Moreover, even when only contacts between children of different classes are taken into account, Table 2 shows that each class has more contacts with the other class of the same grade than with any other class, as also reported in [START_REF] Conlan | Measuring social networks in British primary schools through scientific engagement[END_REF]. Finally, teachers have sparse contacts with one another because they spend most of the time in class with children.

Temporal evolution of the contact patterns

Figure 6 shows the time evolution of the 20-minute aggregated networks for each day of the study. The number of individuals is stable during teaching hours, morning and afternoon breaks, and drops during the lunch time as some children were going back home to have lunch (the school is located in an urban area and many children actually live nearby and can go home for lunch). The average degree of the network displays a more interesting behavior, as it peaks at various moments, each corresponding to a break or the beginning or end of the lunch of a series of classes.

Comparison between day 1 and day 2

A comparison between the characteristics of the overall face-to-face contact patterns in the two days of the deployment is reported in Table 4. Statistical quantities such as the average total number and durations of contacts, the number of different persons contacted, or the contact durations are extremely close across the two days. As shown in Figure 2

Discussion

To our knowledge, this is the first study presenting detailed measures of close (face-toface) proximity interactions between children in a primary school (see however [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF] for the case of a high school). These descriptive results on contact patterns are of interest for modeling the spread of various infectious diseases, and possibly for investigating the role of specific control measures, such as closure of classes, immunization strategies, and so on. Time-resolved contact data were collected over two school days by deploying a wireless sensor network of RFID badges that record close-range (1 to 1.5m) face-to-face proximity between individuals. The present study had a very high participation rate (> 95%). Relying on unobtrusive wearable devices allows the unsupervised detection of contacts during which a communicable disease, in particular a respiratory disease, may be transmitted. This is an important advantage compared to approaches based on questionnaires, especially among the youngest.

Comparison with previous studies

A number of other studies describe or estimate social contact numbers and durations [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF][START_REF] Mikolajczyk | Social contacts of school children and the transmission of respiratory-spread pathogens[END_REF][START_REF] Glass | Social contact networks for the spread of pandemic influenza in children and teenagers[END_REF][START_REF] Zagheni | Using time-use data to parametrize models for the spread of close-contact infectious diseases[END_REF][START_REF] Wallinga | Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents[END_REF][START_REF] Mossong | Social contacts and mixing patterns relevant to the spread of infectious diseases[END_REF][START_REF] Beutels | Social mixing patterns for transmission models of close contact infections: exploring self-evaluation and diarybased data collection through a web-based interface[END_REF][START_REF] Hens | Mining social mixing patterns for infectious disease models based on a two-day population survey in Belgium[END_REF]. Comparison with previous results is clearly important but is made difficult by differences in the definitions of interaction/contact as well as by differences in the measurement techniques. As the present study considers the unsupervised detection of face-to-face proximity, it does not rely on surveys nor on the memories of participants. It is thus expected that larger total number and durations of contacts will be obtained, in comparison with survey-based methods.

Table 5 reports the comparison of the number and duration of contacts between previous studies and the present one. As expected, when all contacts are taken into account, we obtain larger values than the studies cited above, with the exception of [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF]: as the infrastructure described in [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF] considers a broader detection range (3 meters proximity) than in the present case, it is not surprising that our study detects less numerous and shorter contacts. We report that each child has on average 323 contacts lasting 33 seconds per day with other children, corresponding to contacts with an average of 47 distinct other children, for an average daily total interaction time of 176 minutes. The present study gives however access to much more detailed information such as the duration of each contact and the cumulative duration of the contacts between two

individuals. An important result concerns the fact that most contacts are short but that the distribution of contact durations is very heterogeneous, with a non-negligible fraction of longlasting contacts. Strongly heterogeneous distributions of contact durations such as the one displayed in Figure 3 have been observed in other settings, including conferences involving only adults [START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF][START_REF] Alani | Live Social Semantics, 8th International Semantic Web Conference ISWC[END_REF] and a US high school [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF]. Observing similar patterns among young children was however not a priori expected.

To allow a more informed comparison between studies based on different methodologies, we compute for each child or for each pair of individuals the number and total duration of contacts lasting longer than a given threshold. The results are summarized in Table 6. For instance, when restricting to cumulated contact durations of at least one minute, the number of different children with whom a child has interacted drops to 21, and the corresponding total interaction time drop to 163 minutes. Moreover, numbers close to those reported by Glass et al.

[9], Zagheni et al. [START_REF] Zagheni | Using time-use data to parametrize models for the spread of close-contact infectious diseases[END_REF] and Del Valle et al. [START_REF] Valle | Mixing patterns between age groups in social networks[END_REF] are obtained when one takes into account only pairs of children having interacted for at least 10-12 minutes per day. Overall, our results are therefore quantitatively different from other studies, as can be expected from the strong methodological differences, but become compatible with previous studies when applying filtering procedures which retain only the longest contacts.

Public health implications in the field of infectious diseases

Our results show that children mix preferentially with children within their age group.

This effect, known as age homophily, is largely due to the fact that children study together and have the same schedule, and represents a general feature studied in various contexts by sociologists [START_REF] Mcpherson | Birds of a Feather: Homophily in Social Networks}[END_REF]. As a result, the contact matrices display a hierarchical block-diagonal structure once their entries (classes) are sorted according to the seniority of children, as visible in Figures 4 and5 (Tables 2 and3): in addition to the strong diagonal entries (which correspond to the clusters of Figure 9, i.e., to contacts within classes) one can see blocks corresponding to the same grade (e.g., 1A/1B, 2A/2B, 3A/3B, etc.) and two larger blocks that span the junior grades (1 to 3) and the senior grades (4 to 5). In addition, the study of the daily evolution of the contact patterns allows to detect periods of higher contact activity, corresponding to breaks or lunch time, as well as moments when the contact weights (i.e., the cumulative time in contact) increase the most.

These results may help to advise public decision-makers on interventions aimed at containing or mitigating the propagation of communicable diseases at the level of schools, in particular in case of an epidemic or a pandemic. School closure has been proposed as an effective physical intervention to reduce transmission of respiratory pathogens, especially influenza [START_REF] Cauchemez | Closure of schools during an influenza pandemic[END_REF].

However, it is not well understood how the benefit of closing entire schools, in terms of reducing cases, morbidity and mortality, compares to the economic costs of such interventions. In addition, the effectiveness of school closure depends on the effectiveness of other measures such as vaccination or antiviral drugs. Our results could be of interest in this context, especially if combined with other sources of information on the contact patterns of children [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF][START_REF] Mikolajczyk | Social contacts of school children and the transmission of respiratory-spread pathogens[END_REF][START_REF] Conlan | Measuring social networks in British primary schools through scientific engagement[END_REF]. The fact that a child spends three times more time in contact with classmates than with children of other classes suggests for instance that closing selected classes instead of the whole school could be a viable alternative. Additional intermediate steps between class and school closures could be devised through the analysis of aggregated contact networks, such as the one depicted in Figure 8, and exposure matrices such as the ones of Tables 2 and3: classes most strongly linked to the class of the first detected case could for instance be closed in order to reduce the risk of propagation to the remaining classes. It would be interesting to assess by means of numerical simulation whether the closure of a single class or of a group of classes could efficiently mitigate the propagation of a disease at the school level. Finally, as highlighted by Figures 8 and10 and by the movie in the Supplementary Information, preventive measures such as shifts of the class schedules could substantially reduce contacts between classes, which could be particularly relevant for preventing transmission events from asymptomatic cases.

Informing mathematical models

The development of mathematical models that aim at describing the spread of the infection and its prevention and control is hindered by the lack of information on the contact patterns between individuals. Epidemiological models of disease transmissions in structured populations depend heavily on the knowledge of the amount and duration of contacts between individuals of different age groups. To reduce this knowledge gap, we provide the exposure matrix of Table 2. In this matrix, the cell at row A and column B gives the average number and duration of contacts between one individual of category A with any individual of category B. This may help to refine the young age groups of the contact matrices proposed by Zagheni [START_REF] Zagheni | Using time-use data to parametrize models for the spread of close-contact infectious diseases[END_REF],

Wallinga [START_REF] Wallinga | Using data on social contacts to estimate age-specific transmission parameters for respiratory-spread infectious agents[END_REF] and Mossong [START_REF] Mossong | Social contacts and mixing patterns relevant to the spread of infectious diseases[END_REF]. Given the important role of children and young adolescents in the community spread of respiratory infections, it is important to detail this part of the matrix. In particular, we remark that most contacts occur within a given class, and relatively few contacts occur across classes, effectively limiting the ability for diseases to spread between different classes.

These results highlight important properties of the contact patterns between school children that need to be taken into account when modeling the propagation of diseases and when evaluating control measures. On the one hand, our results tend to indicate that assumptions such as the homogeneity of contact durations, or a homogeneous mixing between classes, may yield misleading results. On the other hand, Figures 7 and8 (together with the movie in the Supplementary Information) show that homogeneous mixing within each class seems to be a good approximation. Of course, these considerations should be assessed by numerical simulations comparing the relative outcomes of different assumptions, such as global homogeneous mixing, homogeneous mixing within each class, and a network model.

Caveats

In the following, we discuss some limitations of the present study, and point to strategies for moving forward.

First of all, the deployed infrastructure only measured contacts between children while they were in the school building or in the playground. Badges were not worn during sport activities, which often involve close proximity situations and physical contacts. Moreover, even though the children would not be in school during a school closure, they would still mix with other children and adults in the community and spread the virus through these contacts. It would be interesting to use the data collection infrastructure to combine school data with household data and data on contact patterns during school closure [START_REF] Jackson | School closures and student contact patterns[END_REF][START_REF] Eames | The impact of school holidays on the social mixing patterns of school children[END_REF]. Coupling the dynamical contacts patterns at school and at home would allow to improve our understanding of the role of children as a reservoir in the community spread of infections.

Another potential issue concerns the possibility that children changed their behaviour because they were wearing badges and knew they were participating in a scientific measure.

According to observers familiar with the environment (teachers and staff), however, no significant change could be detected in the children's behavior, and the children seemed to rapidly forget about the badges. In addition, while detailed explanations were given to the parents about the study and the badges, details on the role of the RFID badges (e.g., their detection range)

were not given to the children.

From a public health perspective, it has to be emphasized that the collected data provide information on the mutual proximity of badges (and therefore of the persons wearing the badges), but not on the occurrence of physical contacts. Our measurements may thus be used in the context of, e.g., respiratory-spread pathogens but not for infectious agents transmitted by skin contact. Note however that physical contact can only occur between persons who are already in spatial proximity. Therefore, it would be very interesting to study the fraction of close encounters that result in a physical contact. In the future, the use of devices that can directly sense physical contact (e.g., body-area networks) may be explored.

The short period of time (two days) of data collection also limits the ability to draw conclusions on what happens at longer time scales. Deployment of the sensing infrastructure over much longer timescales is needed in order to confirm the present results.

Finally, the data presented in this study depend on the school schedule and spatial structure, and the generalization of our results to other schools should be carried out with caution.

Some properties are however expected to be rather general, such as the heterogeneity of the contact durations and of the cumulated contact durations, that has been observed in several other settings [START_REF] Salathé | A High-Resolution Human Contact Network for Infectious Disease Transmission[END_REF][START_REF] Hui | Pocket switched networks and human mobility in conference environments[END_REF][START_REF] Cattuto | Dynamics of person-to-person interactions from distributed RFID sensor networks[END_REF][START_REF] Isella | What's in a crowd? Analysis of face-to-face behavioral networks[END_REF][START_REF] Isella | Close encounters in a pediatric ward: measuring face-to-face proximity and mixing patterns with wearable sensors[END_REF] and seems to be ubiquitous in human interactions.

Perspectives

Further research will use the gathered data to simulate the transmission of infectious agents (e.g. respiratory or gastro-enteric viruses) inside a school, to evaluate the role of the index case, and to assess the impact of various containment measures (e.g. class closure, homogeneous partial vaccination vs. vaccination of whole classes, at fixed coverage, etc). Further deployments in other schools with different schedules, other countries, and possibly for longer periods, will also be very useful to cross-validate our findings. ). 88% of the contacts last less than 1 minute, but more than 0.2% last more than 5 minutes. For the cumulated durations, 64 % of the total duration of contacts between two individuals during one day last less than 2 minutes, but 9% last more than 10 minutes and 0.38% more than 1 hour. The small symbols correspond to the original distributions, and the large symbols to the log-binned distributions.

  indicates a high-variance distribution.

--

  the frequency ij n gives the number of times that a contact event between i and j was recorded the cumulative duration ij w is the sum of the durations of the ij n contacts between i and j (

1 )-- 1 ( 1 (

 111 , and the average number of contacts of a child of class A with children of class B, the total time spent in contact between children of classes A and B the average contact time of a child of class A with children of class B, the total number of contacts of children of class A, ∑ B AB n , and the average number of contacts of a child of class A, where the sums on B include the case B=A); the total contact time of children of class A ∑ B AB w , and the average time spent in contact by a child of class A, where the sums on B include the case B=A). The quantities AB n and AB w define symmetric contact matrices at the class level, while the quantities symmetric matrices taking into account the different class sizes

Figure 1

 1 Figure 1 reports the total number and cumulated duration of contacts involving children of a specific class or teachers. Figure 2 displays boxplots of the distributions of the individual contact numbers and cumulated durations, for each class and for each day.Figure 1 shows that

  Figure1shows that the total number and duration of contacts involving teachers are smaller compared to those involving children, but Figure2indicates that this is mostly a consequence of the number of teachers being smaller than the number of children in a class. Both the number and the duration of contacts show a limited degree of heterogeneity across classes as well as across days. This is probably due to different class schedules (e.g., a class being absent during half a day because of sport activities) or to different school activities (e.g., group vs individual activities).Each individual, on average, was in contact with 50 distinct individuals ( ) in day 2. The corresponding distributions are shown as boxplots in Figure2C, for each class and for each day.

Figure 3 Figures 4

 34 Figure 3 displays the distribution of the total amount of time w ij that pairs of individuals i and j spent in contact during one day. Once again, most cumulated durations are short, but their distribution is broad: 64% of the pairs of individuals have interacted less than 2 minutes on the same day, 9% have spent more than 10 minutes together and 0.38% more than 1 hour. The average amount of time spent by two persons in face-to-face proximity during one day is 207 seconds (3min 27s) for day 1, and 236 seconds (3min 56s) for day 2, with squared coefficients of variation of the distributions of respectively 5.44 and 4.65. When considering single individuals, the distribution of total time spent by an individual in face-to-face proximity with other individuals is more homogeneous, with an average of 10340 seconds (2h 52mn) for day 1 and 11000 seconds (3h 03mn) for day 2, with 25 . 0 2 = CV for day

Figure 7

 7 Figure 7 displays the time evolution of the quantities i k , in i k and out i k

Figure 8 Figure 9

 89 Figure 8 gives more insight into the evolution of the contacts of the children by taking into account the cumulated time spent in contact. It shows the time evolution of i s , in i s and

  Figure 10 displays the trajectories followed by children as they move across the classes and public spaces of the school. It shows how each class moves from its classroom to the courtyard and then comes back at various times. During the lunch break, some children go first to the cafeteria and then to the courtyard, encountering children who are moving in the opposite direction. It is apparent how these trajectories, dictated by the school schedule, strongly contribute to shaping the mixing patterns between classes and grades.

Figure 1 .

 1 Figure 1. Total cumulated duration (A, in hours) and number of contacts (B) involving individuals of each class, for each day.

Figure 2 .

 2 Figure 2. Boxplots of the distributions of the cumulated duration (A, in minutes) and number (B) of contacts involving an individual, for each class and for each day. Panel C gives the distributions of the number of distinct individuals with whom an individual of each class has had at least one contact. In each boxplot, the horizontal line gives the median, the box extremities are the 25 th and 75 th percentiles, and the whiskers correspond to the 5 th and 95 th percentiles.

Figure 3 .

 3 Figure 3. Log-log plot of the distribution of the contact durations and of the cumulated duration of all the contacts two individuals i and j have over a day ( ij w ). 88% of the contacts last less

Figure 4 .

 4 Figure 4. Grayscale-coded contact matrix between classes. The matrix entry for row A and column B gives the number of contacts ( AB n ) measured between individuals of classes A and B over the two days of data collection. A logarithmic grayscale is used to compress the dynamic range of the matrix entries and enhance the off-diagonal hierarchical structure.

Figure 5 .

 5 Figure 5. Grayscale-coded contact matrix between classes. The matrix entry for row A and column B gives the cumulated duration ( AB w , in minutes) of contacts measured between individuals of classes A and B over the two days of data gathering. A logarithmic grayscale is used to compress the dynamic range of the matrix entries and enhance the off-diagonal hierarchical structure.

Figure 6 .

 6 Figure 6. Degree of individuals in the contact networks aggregated over sliding time windows of 20 minutes during the first day (left) and the second day (right) of data collection. The median value is represented with a black line, the 95 % confidence interval is shown in gray and the number of individuals over which the statistics are calculated is shown in red dashes. Breaks and

Figure 7 .

 7 Figure 7. Time evolution of the average number of distinct children with whom a child has been in contact during the study. The average total number is displayed in black, the average number of children of the same class in red, and the average number of children of other classes in blue.

Figure 8 .

 8 Figure 8. Time evolution of the average cumulated time spent by a child in contact with other children during the study. The average total time is displayed in black, the average time spent with children of the same class in red, and the average time spent with children of other classes in blue.

Figure 9 .

 9 Figure 9. Network of contacts aggregated over the first day. Edges between individuals having interacted less than 2 minutes have been removed, thus keeping only the strongest links. The width of links corresponds to the cumulative duration of contacts, and nodes with higher number of edges have larger size. Colors correspond to classes, teachers are shown in grey. Figure created using the Gephi software, http://www.gephi.org.

Figure 10 .

 10 Figure 10. Approximate spatiotemporal trajectories of some classes. Each row corresponds to a particular place in the school (classroom, canteen, courtyard) where a RFID reader was situated, and each colored line corresponds to the spatio-temporal trajectory of the children of a class (only 5 classes are shown for clarity). Line widths correspond to the number of children whose approximate position correspond to the row area. A line can become thinner if children leave the school (for instance during the lunch break, to have lunch at home) or divide itself into two thinner lines if two groups of children of the same class follow distinct paths in the school The trajectories highlight how mixing between classes, shown by the fact that the colored lines overlap, occurs during the breaks and is strongly constrained by the school schedule.

Table 1 .

 1 Rate of participation for the school classes.

	Grade	Class name	Numbers of children or teachers	Number of participating children or teachers Day 1 Day 2	Participating rate (%) Day 1 Day 2
	1	Class 1A Class 1B	24 25	22 25	23 25	91.7 100	95.8 100
	2	Class 2A Class 2B	25 26	22 25	23 26	88.0 96.2	92.0 100
	3	Class 3A Class 3B	24 22	23 21	23 21	95.8 95.5	95.8 95.5
	4	Class 4A Class 4B	23 24	21 22	21 22	91.3 91.7	91.3 91.7
	5	Class 5A Class 5B	24 24	22 23	21 23	91.7 95.8	87.5 95.8
	-	Teachers	10	10	10	100	100

Table 3 .

 3 Exposure matrix between grades. The cell of row A and column B of the matrix gives the average number (and the duration in minutes, between parenthesis) of contacts involving an individual of grade A with any individual of grade B, per day.

		1 st grade	2 d grade	3 d grade	4 th grade	5 th grade	Teachers
	1 st grade	322.6 (177.7)	24.8 (13.3)	13.9 (6.2)	10.0 (4.7)	10.4 (4.4)	13.0 (5.5)
	2 d grade	24.6 (13.2)	274.8 (162.7)	21.7 (11.4)	3.0 (1.2)	2.1 (0.7)	11.4 (6.0)
	3 d grade	15.0 (6.7)	23.9 (12.5)	307.7 (167.8)	7.7 (3.5)	6.4 (2.4)	7.9 (4.0)
	4 th grade	11.1 (5.2)	3.3 (1.4)	7.9 (3.6)	189.9 (103.7)	18.2 (9.2)	6.4 (2.9)
	5 th grade	11.3 (4.8)	2.3 (0.8)	6.3 (2.4)	17.5 (8.9)	269.3 (148.6)	7.3 (3.6)
	Teachers 61.7 (26.1)	54.8 (29.0)	34.8 (17.5)	27.5 (12.4)	32.4 (16.2)	10.5 (6.8)

Table 4 .

 4 Comparison of some characteristics of the networks of day 1 and 2.

		Day 1	Day 2
	Number of individuals	236	238
	Average number
	of contacts of an	317 (0.22)	338 (0.27)
	individual (CV 2 )
	Average total
	time in contact of an individual, in	172 (0.25)	183 (0.33)
	minutes (CV 2 )
	Average number
	of distinct persons	50 (0.14)	46.5 (0.18)
	contacted (CV 2 )
	Average
	cumulated time
	spent in contact	207 (5.4)	236 (4.7)
	by two persons,
	in seconds (CV 2 )
	Average duration
	of a contact , in	32.6 (1.2)	32.6 (1.1)
	seconds(CV 2 )
	Average
	clustering	0.5	0.56
	coefficient

Table 6 .

 6 A. Average number of children with whom a child was in contact, computed over one day, and average total time spent daily in contact with other children. In both cases, a filtering that retains only the contacts with duration at least equal to T is applied (T=0 or 20s corresponds to taking all contacts into account, given the available 20s time resolution). B. Average number of children with whom a child was in contact, computed over one day, and average total time spent daily in contact with other children. For both quantities, a filtering procedure is applied that retains only the links between children who have spent an amount of time at least equal to W in face-to-face proximity (W=0 or 20s corresponds to taking all links into account).

	A. Filtering procedure: only contacts	Average daily number of distinct	Average daily cumulated duration of
	of duration at least T	other children in contact	contacts with other children, in
			minutes
	T=0	47.4	176
	T=40s	20.8	100
	T=1 mn	11.8	65
	T=2 mn	4..1	28
	T=3 mn	2.2	19
	B. Filtering procedure: only	Average daily number of distinct	Average daily cumulated duration of
	cumulated contacts at least W	other children in contact	contacts with other children, in
			minutes
	W=0	47.4	176
	W=1 mn	21.4	163
	W=2 mn	15.2	153
	W=5 mn	8.1	129
	W=7 mn	6.1	117
	W=10 mn	4.3	102
	W=12 mn	3.5	93
	W=15 mn	2.7	81

Acknowledgments

We warmly thank Bitmanufaktur, the Openbeacon project and Truelite for their technical support. We are particularly grateful to all children, their parents and the school staff who volunteered to participate in the data collection.

Setting

Contact definition Results

Mikolajczyk et

al. [START_REF] Mikolajczyk | Social contacts of school children and the transmission of respiratory-spread pathogens[END_REF] Survey in a primary school; 6-10 yearsold children.

A person with whom the child spoke or played with in a day 25.1 contacts per day per child -Wallinga et al.


Jérôme Schmid
email: schmid@miralab.unige.ch

Anders Sandholm
email: anders.sandholm@epfl.com

François Chung
email: francois.chung@sophia.inria.fr

Daniel Thalmann
email: daniel.thalmann@epfl.com

Hervé Delingette
email: herve.delingette@sophia.inria.fr

Nadia Magnenat-Thalmann

Musculoskeletal simulation model generation from MRI datasets and motion capture data

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

Introduction

Today neuromuscular simulation has proved to be a valuable tool in understanding human movements, from gait research [START_REF] Piazza | The influence of muscles on knee flexion during the swing phase of gait[END_REF], diagnostic and treatment of patients to teaching/training new physicians. During the last years, several simulation platforms have been developed, both commercial [START_REF] Damsgaard | Analysis of musculoskeletal systems in the anybody modeling system[END_REF][START_REF] Delp | An interactive graphicsbased model of the lower extremity to study orthopaedic surgical procedures[END_REF] and open source based [START_REF] Delp | Opensim: Opensource software to create and analyze dynamic simulations of movement[END_REF], along with several models with different level of details and complexity [START_REF] Zhang | Dynamic and static control of the human knee joint in abductionadduction[END_REF][START_REF] Holzbaur | A model of the upper extremity for simulating musculoskeletal surgery and analyzing neuromuscular control[END_REF].

The use of computer model/simulation has also allowed researchers and physicians to test what if scenarios, and even simulate different treatments before the physical treatment or surgery take place. A person gait cycle can be affected by many factors such as bone deformations [18] or static muscle contractors [START_REF] Tardieu | Mechanical evaluation and conservative correction of limb joint contractures[END_REF]. To simulate such patient with gait troubles or deformed skeleton, a full inspection of their anatomy should be first conducted. The most common neuromuscular models which are used today are all based on healthy, average-size adult male specimens [START_REF] Horsman | Morphological muscle and joint parameters for musculoskeletal modelling of the lower extremity[END_REF]. Before the simulation of a person anatomy, the generic neuromuscular model has to be scaled to match the anatomy of the subject. The most two common ways are the isotropic scaling, where the whole model is scaled with the same scaling factor and the anisotropic scaling where each model part is scaled with an individual scaling factor. To determine these scaling factors, an initial motion capture is performed, where the subject stands still wearing 3D position markers placed at pre-defined anatomical positions. From these 3D markers, the length and position of each body segment can be computed to estimate the appropriate scaling factors to be applied on the generic musculoskeletal model. This scaling approach has several limitations and drawbacks. Firstly, the surface-based recordings cannot fully account for the large anatomical differences that exist among individuals. Secondly, they do not tackle the issue of skin displacement with respect to the underlying bones. A slight misplacement of a marker can generate large differences in the simulated model. In patient with cerebral palsy, it has been shown that the length of moment arms can significantly differ between a scaled and a subject-specific model [START_REF] Scheysa | Personalized mrbased musculoskeletal models compared to rescaled generic models in the presence of increased femoral anteversion: effect on hip moment arm lengths[END_REF].

The scaling is not the only issue; the complexity of the musculoskeletal system should be also addressed. Various studies tried to set up anatomical models at various levels of complexity. Scheepers et al. [START_REF] Scheepers | Anatomy-based modeling of the human musculature[END_REF] and Aubel and Thalmann [START_REF] Aubel | Interactive modeling of the human musculature[END_REF] described approaches based on anatomical concepts and constraints but the methods used unrealistic simplifications. One of the major limitations was that both studies were not subject-specific. More advanced models were presented like the works of Teran et al. [START_REF] Teran | Creating and simulating skeletal muscle from the visible human data set[END_REF] or Blemker and Delp [START_REF] Blemker | Three-dimensional representation of complex muscle architectures and geometries[END_REF] creating more complex anatomical models (e.g., fiber directions) to increase simulation accuracy. Although these studies relied on medical data sets offering the individualization of the modeling, they still require manual operations and corrections.

In this paper, a new framework that combines data from motion capture and segmented MRI images is introduced to create a subject musculoskeletal model and therefore remove the initial scaling step. The anatomical modeling pipeline uses first a stage of interactive modeling to create what is referred to as a generic anatomical model. This modeling needs to be performed only once and it considers anatomical constraints. Then the generic anatomical model is used in an automatic procedure that can be applied on any dataset to perform the individualization. The proposed approach can thus relieve the user from tedious interactive modelings. Finally, the anatomical model is combined with additional information produced from a reference motion capture frame to generate the complete musculoskeletal model in a semi-automatic fashion. If needed, a user interface is provided to alter/refine the full model. This model can be then used to simulate various motions.

Musculoskeletal modeling pipeline

When generating a musculoskeletal model, several sources of information and data have to be combined in a single simulation setup. The overall workflow, depicted in Fig. 1, is divided up into two phases, a and b. In phase a, an individualized musculoskeletal model describing the subject is created, this model can be subsequently used in phase b where a neuromuscular simulation is carried out to analyze a desired motion. In the first phase, a data acquisition process takes place. The lower limbs of the subject are scanned in a MRI device to produce datasets that are later on segmented (see corresponding section) in various anatomical components (e.g., muscles, bones). Motion capture sessions are also performed to collect information during a standing reference motion capture and one or several experiment motion captures. The later motion capture is coupled with synchronized EMG acquisitions and is used to simulate the motion that we want to study. By combining the 3D models produced from the segmentation with the reference motion capture data, the subject-specific musculoskeletal model is generated. This model is then used in the second phase, along with the experiment motion capture data to simulate and analyze a specific motion.

In developing this framework, the OpenSim platform was chosen as the neuromuscular simulation engine, mainly because it is based on an open-source license and it supports several tools for model editing, fast simulation and advanced visualization features. OpenSim also supports a large collection of simulation possibilities such as inverse kinematics, muscular control and forward dynamics, whose output parameters include joint forces, muscle-induced acceleration and muscle power.

Data acquisition

Markers setup for MRI and motion capture

To create the musculoskeletal model, two marker sets are utilized. One set of skin markers is placed on muscle tissue while the other set of reference markers is placed on bone features, see Table 1 and right side of Fig. 1. The reference markers are located on easy-to-find bone features and are used both during the MRI acquisition and the motion capture. During motion capture, reflection markers are used, while markers filled with contrast agent are used during the MRI acquisition. The reference markers are used to put in correspondence the segmented and the motion capture data. Due to the fact that the reference markers are placed by hand at two different acquisitions, a small difference in placement can occur. However, this error is smaller than the error produced by skin movements during acquisition.

MRI acquisition

In close collaboration with radiologists, an adequate protocol for the imaging of soft and bony tissues was defined: Axial 2D T1 Turbo Spin Echo, TR/TE = 578/18ms, FOV/FA = 40cm/90°, matrix/resolution = 512x512/0.78x0.78mm. To acquire a full lower limb, two different sessions took place and during each of them three acquisitions were performed. Each acquisition was based on the same protocol, but with varying slice thickness (e.g., higher thickness for long bones) to speed up the acquisition process. First session consisted in hip (thickness: 2mm), thigh (10mm) and knee (4mm); second session in knee (2mm), leg (10mm) and foot (4mm). All acquisitions of a session were merged in a unique volume, and the two session volumes were registered together (thanks to an appropriate knee overlap). After cropping and resampling operations, two datasets have been created with the following size (resolution): 202 x 398 x 595 (0.78 x 0.78 x 2.46mm) and 151 x 213 x 582 (0.68 x 0.68 x 2.5mm) for the first and second datasets respectively. All acquisitions were performed at the University Hospital of Geneva on a 1.5T MRI device (Philips Medical Systems).

Motion capture acquisition

For the recording of the reference and experiment motion capture data, an optical motion capture system of eight ProReflex Qualisys (www.qualisys.com) video cameras was used. Prior to the session, reflective skin markers were fitted to the lower limbs according to Table 1. Subjects were also equipped with electromyography (EMG) sensors from Noraxon (www.noraxon.com), to capture electrical activity in eight large lower limb muscles, see Table 2. During motion acquisition, two AMTI (www.amti.biz) force platforms were used to capture ground reaction force. Marker trajectories, EMG and force plate data were acquired, labeled and exported by using Qualisys Track Manager to C3D files.

Segmentation

In order to model the various structures of interest with a high fidelity, a segmentation procedure is applied on medical datasets of the subject. The usage of manual segmentation is twofold. Firstly, it serves as a basis to automatic methods by providing (i) some prior knowledge (e.g., internal forces in deformable models that enforce the consistency with an a priori of a shape [8,[START_REF] Gilles | Anatomical modelling of the musculoskeletal system from mri[END_REF]) or (ii) training material (e.g., methods based on Principal Component Analysis (PCA) [START_REF] Cootes | The use of active shape models for locating structures in medical images[END_REF][START_REF] Schmid | MRI bone segmentation using deformable models and shape priors[END_REF] or classification methods like K-means [START_REF] Hartigan | A K-means clustering algorithm[END_REF] and Fuzzy C-means [START_REF] Bezdek | Fuzzy mathematics in pattern classification[END_REF]). In general, more than one model is needed to produce satisfactory results. Secondly, it provides a ground truth for experiments. Segmentations generated by an automatic method can be thus compared with those obtained manually. In our case, manual segmentation has been performed by a medical expert, which allows us to validate it. In the following, our methodology for manual segmentation will be presented then a semi-automatic segmentation method will be depicted.

Generic anatomical model construction

Manual delineation

The manual delineation of the bones and muscles has been performed by a medical expert (see Fig. 2a). For each structure, a binary image was produced and then processed by the Marching Cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF] to generate 2-simplex meshes [8] representing the structures of interest. However, the generated meshes surface appeared not to be smooth enough. This is mostly explained by the dataset resolution (especially the large interslice distance) and the ubiquitous error made during a manual segmentation, which remains a subjective task. Ideally, the resolution should be as high as possible to provide detailed images and to avoid big jumps between consecutive MRI slices. This is not always feasible (e.g., device limitations, acquisition time restrictions). Furthermore, a lower resolution is preferred to reduce the memory footprint which results in a speed-up automatic algorithm. Since muscles are known to be anatomically smooth, the objective is to find a solution that makes them appear so. The technique we used consists first in refining the meshes and then in applying internal constraints on them. More precisely, we applied forces on the meshes so that their rigidity is the maximum. For that, we used the simplex angle continuity or 2 C constraint [8] used for shape recovery and which is defined as an average of the simplex angles at neighborhood vertices. Since we want to smooth our meshes in a large scale, we defined a great neighborhood size which intuitively corresponds to the notion of high rigidity for the mesh. The result was that most pits and bumps present on the meshes surface disappeared throughout the deformation process (see Fig. 2c). The smoothing based on 2 C constraint has also the advantage to reduce the shrinking effect obtained with standard Laplacian smoothing.

Tendons

Although the muscles are difficult to delineate manually, the tendons are easier to detect in the MRI. Indeed, they have a uniform appearance in the images (tendons appear as dark structures with respect to muscle tissue). An automatic method based on maximum intensity ridge tracing [START_REF] Aylward | Initialization, noise, singularities, and scale in height ridge traversal for tubular object centerline extraction[END_REF]20] was first used. Such method relies on the assumption that tendons are roughly tubular structures, which is especially the case for the leg tendons. These tracing methods have proved to be robust to noise. They provide also a simple way to model a structure as a series of centerline points (a position and a radius). Furthermore, such representation could be reused in another automatic method. However, this method may not always segment correctly or entirely the tendons because in some regions their intensity is corrupted by too strong artifacts (e.g., noise, adjacent structures, partial volume effect). As a result, tendons are difficult to follow through the slices. Moreover, some tendons are so close (e.g., fibularis longus and brevis tendons) that it is not possible to distinguish them without some prior knowledge. To cope with these problems, missing centerline points were manually placed (Fig. 3a) on high curvature points of tendons using an editing tool, such as Car-dioViz3D [START_REF] Toussaint | An Integrated Platform for Dynamic Cardiac Simulation and Image Processing: Application to Personalised Tetral-ogy of Fallot Simulation[END_REF]. Eventually, when the centerlines are complete (thanks to the automatic or manual modeling), simplex surfaces of tendons are automatically created. Fig. 3b depicts some of the results for the foot tendons complex modeling.

Attachments

Anatomically, an attachment is defined as the linking region between two structures. For instance, tendons are usually attached to bones. The tendon attachments are thus defined as the tendon tissues sharing a common region with the bone tissues. Comparatively, it means that the meshes modeling these structures should be stuck together which is not the case a priori. Indeed, as accurate as the segmentation may be, the different generated meshes are not guaranteed to be attached together, especially after the smoothing procedure. Our solution is to manually place splines to define attachment regions as described in [START_REF] Gilles | Anatomical and kinematical modelling of the musculoskeletal system from mri[END_REF]. This provides a simple but efficient way to model attachment regions. Meshes are then deformed until they are stuck.

Interpenetrations removal

The smoothing and attachment procedures previously described can create non realistic interpenetrations between meshes. In order to remove them, a geometrical post-processing is applied. Let's consider the surfaces 1 S and 2 S in Fig. 4 that are interpenetrating each other. The aim is to move the colliding points of each surface, located in the penetration volume, to reach a non interpenetrating state. S in e . This point e must be inside 1 S to be valid. The new position m of the point s is then chosen as: S . This can be useful to resolve penetrations between a bone (very stiff material) and a muscle. The presented technique was inspired from [START_REF] Keiser | Contact handling for deformable point-based objects[END_REF], with the distinction that point e is not defined as the projection of s on 2 S . Indeed, by using this projection, the point m does usually not correspond to the new position of s in the non colliding state. By using instead the normal direction n , the point m will provide a better approximation. As a result, this technique avoids the usage of a resampling scheme as described in [START_REF] Keiser | Contact handling for deformable point-based objects[END_REF], points just need to be moved to their new positions. A perfect contact between the surfaces is nevertheless not ensured but the results are satisfactory when meshes resolution is not too low. Finally, implicit surfaces [START_REF] Alexa | On normals and projection operators for surfaces defined by pointsets[END_REF] are not used since the hypothesis of vicinity to surfaces may be invalid (i.e. in case of big interpenetrations). To speed up the process, precomputed signed distance maps can be used as well as the efficient golden section search technique [START_REF] Press | Numerical Recipes in C[END_REF] to minimize f .

[ ] { } ( ) ⎩ ⎨ ⎧ - = ∈ = 2 2 2 1 1 ) ( * ) ( * ) ( minimal is ) ( / , x d l x d l x f m f e s x m (1)

Semi-automatic segmentation

The result of the previous section is a generic anatomical models collection of the various soft (muscles, tendons) and bony structures, from which various levels of details (LODs) are computed (4 for bones and 3 for soft structures). The objective now is to re-exploit these data as prior models in a semi-automatic method. This method could be then applied on any similar dataset coming from other subjects. The proposed method is mostly based on the deformable models method of Gilles et al. [START_REF] Gilles | Anatomical modelling of the musculoskeletal system from mri[END_REF]. Mesh vertices are considered as lumped mass particles subjected to internal/external forces and evolving under the Newtonian law of motion.

Internal forces regulate shape evolution by enforcing constraints [8,[START_REF] Gilles | Anatomical modelling of the musculoskeletal system from mri[END_REF] on the surface regularity (smoothness, curvature) and on the shape configuration. Shape constraints are simply enforced by predefining their simplex parameters. We added other priors [START_REF] Schmid | MRI bone segmentation using deformable models and shape priors[END_REF] based on PCAs of training shapes and on the modeling of local deformations as a Gauss-Markov random process. External forces attract meshes towards anatomical boundaries by using image information (e.g., intensity, gradient) and impose non-penetration constraints. The attachments use mass modification methods to constrain vertices, while collision handling techniques prevent mesh penetrations through penalty forces.

The Newton equation relates forces to particle state. The resulting differential equations system is time discretized and solved by a stable implicit integration scheme [START_REF] Volino | Implicit midpoint integration and adaptive damping for efficient cloth simulation[END_REF]. The different mesh LODs are used in a multilevel framework using a force propagation mechanism. This strategy confers robustness and accuracy to the segmentation method. Conducted experiments [START_REF] Gilles | Anatomical modelling of the musculoskeletal system from mri[END_REF][START_REF] Schmid | MRI bone segmentation using deformable models and shape priors[END_REF] reported an average ac-curacy of 1.5mm for bone and muscle segmentation. The described segmentation procedure is equivalent to a Model to Image Registration. It is semi-automatic because generic anatomical models are initialized manually by placing landmarks. The power of using and registering generic shapes is that exact geometric correspondences are obtained (i.e., morphological features have the same vertex indexes across individual models, like attachments for instance).

Anatomical modeling validation

The proposed methods for correction (smoothing and interpenetrations removal) and the semi-automatic algorithms are still prone to errors. Indeed, they cannot always guarantee that the segmentation will perfectly delineate the structures of interest. To correct the most significant errors, a medical expert can define manually some constraints by placing points (e.g., points on the organ boundaries) in the MRI datasets. Forces are then applied on the meshes to account the expert constraints [START_REF] Gilles | Anatomical and kinematical modelling of the musculoskeletal system from mri[END_REF]. Similarly, forces based on image information (e.g., gradient) can be used to perform some local corrections. This approach gives nice results but as a general remark, one should avoid to place too many constraints or be tempted to blindly trust the image information. Indeed, the placement of constraints can present the same pitfalls as a manual delimitation: it is difficult to assess the organs transition between slices and ensure meshes smoothness. Consequently, it is essential to find a trade-off between reconstructed models quality and exact segmentation, since images cannot be totally trusted. Most importantly, the purpose of the segmentation is to provide models that must be reusable in our pipeline or in other studies.

Musculoskeletal model generation

The generation of a subject-specific musculoskeletal model can be divided into two parts. First, the anatomical segmented 3D objects have to be aligned with the reference motion capture data and secondly, the segmented objects have to be extended with additional information to create a complete musculoskeletal model.

Data alignment

The first step in generating a subject-specific musculoskeletal model is to align the segmented data with the motion capture reference frame, as depicted in the first step of the workflow (see Fig. 1). Doing this, the correct placement of the mo-tion capture skin markers can be included into the model. This alignment is also performed to prepare the anatomical 3D objects for visualization purpose.

The alignment starts by reading the reference markers from the imported C3D file along with the reference markers from the segmented MRI images. For each bone segment (pelvis, femur and tibia), a rigid transformation matrix is calculated so that each segment can be correctly aligned with the motion capture data. To translate each bone segment, its corresponding markers are used, see Table 3. Currently, the marker extraction in the MRI images is performed manually but an automatic method could be achieved by using standard thresholding and mathematical morphology operators. The combination of various subject-specific data avoids performing any error-prone scaling.

By using each bone segment transformation, the skin markers defined in the reference motion capture file can be similarly expressed in the MRI coordinate space. For example, the skin marker Thigh (see Table 1) is transformed using the Femur matrices (see Table 3).

Model generation

When the initial step of alignment is completed, the model consists of the 3D segmented objects and of motion capture markers. To generate a complete musculoskeletal model, information about muscles and joints has to be included. From the segmented 3D objects, each joint center is calculated. This joint center is used as the center point for the moment force generated by the muscles. Then for each joint, a kinematic function is defined to describe the kinematic behavior of each joint. By using this function and a given joint angle, the relative attachment point of the joints can be determined. Currently, this data is derived from Delp et al. [START_REF] Delp | An interactive graphicsbased model of the lower extremity to study orthopaedic surgical procedures[END_REF]. To calculate the force that each muscle and its tendons can generate, each muscle is expressed into a set of one or more muscle-tendon actuators. These muscle-tendon actuators stretch between the muscles insertion points either as a straight line or via a wrapping point. At a wrapping point, the actuator changes direction to better describe the shape of the muscle as shown in Fig. 5. Wrapping points are also used to prevent muscle-tendon actuators to collide and enter a bone segment. If such a collision occurs, the wrapping becomes a fixed point and the actuator "wraps" around the bone, see Fig. 5b and5c. The muscle-tendon actuator (Fig. 5a) is connected to each joint by using the tendon data from the segmentation process (see section "Attachments"), and the framework supports manual correction of both insertion and wrapping points. Currently, the description for all muscle-tendons and wrapping points are derived from Delp et al. [START_REF] Delp | An interactive graphicsbased model of the lower extremity to study orthopaedic surgical procedures[END_REF] and OpenSim [START_REF] Delp | Opensim: Opensource software to create and analyze dynamic simulations of movement[END_REF]. Eventually, the model is ready to be exported as a set of XML files in the OpenSim format. When a subject-specific simulation is carried out, the generated musculoskeletal model is loaded along with the experiment motion capture data. The model will now match the motion capture data, so no initial scaling has to be carried out. The only initial step before starting a simulation is to align each joint angle in the model with the corresponding joint angle in the first frame of the experiment motion capture. This is done using an inverse kinematic calculation which is performed by the OpenSim framework. The musculoskeletal model should now match the first frame of the motion capture data and the desired simulation can be executed.

Results and discussion

Anatomical modeling results

To illustrate the anatomical modeling, we applied our methodology on datasets produced from a female and healthy subject (age: 24, height: 1m68, weight: 58kg) acquisition. The institutional medical-ethical committee approved the study and the subject gave written informed consent. The manual modeling used to construct the generic anatomical model is applied on the leg area (knee to foot) while the 16 Jérôme Schmid , Anders Sandholm , François Chung , Daniel Thalmann , Hervé Delingette and Nadia Magnenat-Thalmann automatic method is performed on the hip-thigh-knee area. The automatic method relies on generic models of the thigh previously generated from a similar manual interactive method. Results are visible in Fig. 6a to 6d.

Musculoskeletal modeling results

Effort has been put in devising a framework to set up a musculoskeletal model from subject-specific data. Currently, this workflow has only been used on one subject, some results are visible in Fig. 6e. Therefore, no quantitative comparison between scaled and subject specific simulation is available. However, the initial step of scaling has been avoided and this should therefore remove a source of error in the musculoskeletal simulation. The proposed approach also supports the generation of musculoskeletal models from subject with no conventional bone structures (e.g., pathological bones), where scaling would not have been possible since the underlying bone structure can be hard to determine from the surface markers/measurements.

Conclusion and future work

The workflow outlined in this paper allows the user to create a subject-specific musculoskeletal model based on MRI and motion capture acquisitions. The proc-ess is described as a semi-interactive approach where the segmentation is performed using a generic anatomical model. The segmented data is then extended with joint and muscular information so that a full musculoskeletal model is created. Currently, joint kinematics and tendon-muscular actuator information are derived from recognized literature. When the specific musculoskeletal model of the subject is generated, it can be used with any subject motion capture data to produce a desired simulation. The versatility of the framework is one of its qualities as it is not anymore necessary to rely on a standardized healthy male musculoskeletal model. It should be hence possible to simulate patients with disabilities as long as the acquired data reflects their pathology (e.g., organ shapes, abnormal EMG patterns, not conventional motions). Moreover, the framework being coupled with the OpenSim platform, a wide variety of measurements and simulation scenarios are available.

Yet a lot of work remains to fully validate the framework. Especially, more experiments must be carried out to compare the standard approach based on scaling and the proposed one. Accent will be put in studying more carefully the joint kinematics and the muscle-tendon actuators. The proposed workflow will have also to study, from a simulation viewpoint, the differences that exist among individuals such as gender, size and age. The approach will be then used in a medical context where some specific motions on pathological subjects will be examined. We hope that it will give some highlights into the mechanisms of pathologies and their treatments.

Fig. 1 .

 1 Fig. 1. Left: Workflow of subject-specific neuromuscular model generation and simulation. Phase a generates the subject-specific neuromuscular model, in phase b simulation and visualization are performed. Right: Reference (blue) and Skin (red) markers position overview, marker numbers correspond with Table 1. Only the lower limbs markers are shown.

Fig. 2 .

 2 Fig. 2. a) Axial view of the leg with manual delineation in half transparency; b)-c) result of the smoothing procedure on the Flexor Digitorum Longus muscle: (b) mesh surface generated from the manual delineation and (c) mesh surface smoothed using the simplex angle continuity or 2 C constraint.

Fig. 3 .

 3 Fig. 3. Left: manual placement of centerline points for the Extensor Hallucis Longus tendon starting from the eponymous muscle until its attachment above the big toe. Right: results on various foot tendons.

Fig. 4 .

 4 Fig. 4. Top-left: surfaces S1 and S2 are interpenetrating each other. Bottom-left: interpenetrations removal results. Right: results in real cases. It can be seen that muscles are not diffused through the neighbors anymore.

For a point 1 S

 1 S ∈ with inner normal n , the line passing through s and directed by n intersects 2

Fig. 5 .

 5 Fig. 5. a) Visualization of muscle-tendon actuators in Gluteus medius model; b-c) Muscle-tendon actuator representing Gluteus maximus; images b and c show the influence of wrapping points.

Fig. 6 .

 6 Fig. 6. Anatomical and musculoskeletal modeling results: a) Lateral and posterior view of the leg based on manual modeling; b) slice with segmentation contours overlay of the leg; c) slice with segmentation contours overlay of the thigh; d) lateral view of the thigh based on automatic segmentation; e) Musculoskeletal simulation in OpenSim showing bone surface, muscle-tendon actuators and motion capture markers.

Table

  . Only the lower limbs markers are shown.

Table 1 .

 1 Reference and skin markers.

		Placement	Type	Marker Set
	1	Acromioclavicular	Single Marker	Reference Marker
	2	Iliac crest	Single Marker	Reference Marker
	3	ASIS	Single Marker	Skin Marker
	4	Greater trochanter	Single Marker	Reference Marker
	5	Sacrum	Single Marker	Skin Marker
	6	Lateral knee	Single Marker	Reference Marker
	7	Medial knee	Single Marker	Reference Marker
	8	Lateral ankle	Single Marker	Reference Marker
	9	Medial ankle	Single Marker	Reference Marker
	10	5 th metatarsal head	Single Marker	Skin Marker
	11	1 st metatarsal head	Single Marker	Skin Marker
	12	Thigh	Cluster Set	Skin Marker
	13	Shank	Cluster Set	Skin Marker

Table 2 .

 2 Jérôme Schmid , Anders Sandholm , François Chung , Daniel Thalmann , Hervé Delingette and Nadia Magnenat-Thalmann Electromyography electrodes placement.

	Nr	Muscle
	1	Tibialis anterior	with additional ground marker
	2	Gastrocnemius
	3	Vastus medialis
	4	Soleus
	5	Vastus lateralis
	6	Rectus femoris
	7	Biceps femoris
	8	Gluteus maximus

Table 3 .

 3 Markers used to transform each bone segment.

	Bone segment	Reference motion capture marker
	Pelvis	left and right side markers of Iliac crest,
		left and right side markers of Greater tro-
		chanter
	Femur	Greater trochanter, Lateral knee, Medial
		knee
	Tibia	Lateral knee, Medial knee, Lateral ankle
		Medial ankle

Jérôme Schmid , Anders Sandholm , François Chung , Daniel Thalmann , Hervé Delingette and Nadia Magnenat-Thalmann

Jérôme Schmid , Anders Sandholm , François Chung , Daniel Thalmann , Hervé Delingette and Nadia Magnenat-Thalmann

Jérôme Schmid , Anders Sandholm , François Chung , Daniel Thalmann , Hervé Delingette and Nadia Magnenat-Thalmann

Acknowledgments This work is supported by the 3D Anatomical Human project (MRTN-CT-2006-035763) funded by the European Union.