Jiachen Ye

Qitong Hu

Peng Ji
email: pengji@fudan.edu.cn

Marc Barthelemy
email: marc.barthelemy@ipht.fr

The effect of interurban movements on the spatial distribution of population in China

Understanding how interurban movements can modify the spatial distribution of the population is important for transport planning but is also a fundamental ingredient for epidemic modeling. We focus here on vacation trips (for all transportation modes) during the Chinese Lunar New Year and compare the results for 2019 with the ones for 2020 where travel bans were applied for mitigating the spread of a novel coronavirus . We first show that these travel flows are broadly distributed and display both large temporal and spatial fluctuations, making their modeling very difficult. When flows are larger, they appear to be more dispersed over a larger number of origins and destinations, creating de facto hubs that can spread an epidemic at a large scale. These movements quickly induce (in about a week) a very strong population concentration in a small set of cities. We characterize quantitatively the return to the initial distribution by defining a pendular ratio which allows us to show that this dynamics is very slow and even stopped for the 2020 Lunar New Year due to travel restrictions. Travel restrictions obviously limit the spread of the diseases between different cities, but have thus the counter-effect of keeping high concentration in a small set of cities, a priori favoring intra-city spread, unless individual contacts are strongly limited. These results shed some light on how interurban movements modify the national distribution of populations, a crucial ingredient for devising effective control strategies at a national level.

I. INTRODUCTION

In early January 2020, we observed the outbreak of a novel coronavirus (COVID-19) in Wuhan, China, which has been quickly spreading out to the whole country, and more recently to other countries in the world [START_REF] Hernández | declares global emergency as wuhan coronavirus spreads[END_REF]. Infectious diseases spread among humans because of their interactions and movements, and the proximity of this outbreak with the Spring Festival, a period of travel with high traffic loads, provided terrible conditions for the spread of this disease. With an increasing amount of confirmed cases, more attention has been devoted to modeling the spread of COVID-19 from various aspects such as determining the value of the reproductive number [START_REF] Cao | Estimating the effective reproduction number of the 2019-ncov in china[END_REF][START_REF] Tang | Estimation of the transmission risk of the 2019-ncov and its implication for public health interventions[END_REF][START_REF] Riou | Pattern of early human-to-human transmission of wuhan 2019 novel coronavirus (2019-ncov), december 2019 to[END_REF][START_REF] Zhao | Preliminary estimation of the basic reproduction number of novel coronavirus (2019-ncov) in china, from 2019 to 2020: A data-driven analysis in the early phase of the outbreak[END_REF][START_REF] Park | Reconciling early-outbreak preliminary estimates of the basic reproductive number and its uncertainty: a new framework and applications to the novel coronavirus (2019-ncov) outbreak[END_REF][START_REF] Zhang | Evolving epidemiology of novel coronavirus diseases 2019 and possible interruption of local transmission outside hubei province in china: a descriptive and modeling study[END_REF], of the incubation period [START_REF] Jantien | Incubation period of 2019 novel coronavirus (2019-ncov) infections among travellers from wuhan, china, 20-28 january 2020[END_REF][START_REF] Read | Novel coronavirus 2019-ncov: early estimation of epidemiological parameters and epidemic predictions[END_REF][START_REF] Liu | -ncov)[END_REF]. In general, analytical modeling plays of course an important role in the prediction of the spread and allows in particular to test control strategies [START_REF] Heesterbeek | Modeling infectious disease dynamics in the complex landscape of global health[END_REF], which was verified in this case too [12-20, 20, 21, 21, 22]. Particularly important were estimate of probability to export the disease in other countries [START_REF] Chinazzi | Preliminary assessment of the international spreading risk associated with the 2019 novel coronavirus (2019-ncov) outbreak in wuhan city[END_REF][START_REF] Pullano | Novel coronavirus (2019-ncov) early-stage importation risk to europe, january 2020[END_REF][START_REF] Gilbert | Preparedness and vulnerability of african countries against importations of covid-19: a modelling study[END_REF], and how effective were travel restrictions inside China [START_REF] Chinazzi | Preliminary assessment of the international spreading risk associated with the 2019 novel coronavirus (2019-ncov) outbreak in wuhan city[END_REF].

Demographic information and mobility, either under the form of data or given by transportation models (see for example the review [START_REF] Barbosa | Human mobility: Models and applications[END_REF]), are crucial for these transmission models. Mobility can concern either the global scale with movements between countries [START_REF] Joseph T Wu | Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study[END_REF][START_REF] Chinazzi | Preliminary assessment of the international spreading risk associated with the 2019 novel coronavirus (2019-ncov) outbreak in wuhan city[END_REF][START_REF] Lai | Preliminary risk analysis of 2019 novel coronavirus spread within and beyond china[END_REF][START_REF] Chinazzi | The effect of travel restrictions on the spread of the 2019 novel coronavirus (2019-ncov) outbreak[END_REF], or the national scale between cities, or even inside cities [START_REF] Lai | Preliminary risk analysis of 2019 novel coronavirus spread within and beyond china[END_REF][START_REF] Chinazzi | The effect of travel restrictions on the spread of the 2019 novel coronavirus (2019-ncov) outbreak[END_REF][START_REF] Xinwu Qian | Scaling of contact networks for epidemic spreading in urban transit systems[END_REF].

In this study we are focusing on the national level and we won't try to model the spread of the disease. Instead we will focus on the statistical properties of the interurban mobility, and how it affects the spatial distribution of populations. More precisely, we will investigate the statistical properties, of traffic flows between cities during the Chinese Spring Festival in 2020 and by comparing with data for 2019, which are the most salient differences induced by travel restrictions. This knowledge will help us to understand the effect of travel restrictions and their impact on epidemic spread, and more generally to guide us for modeling mobility at this scale, a crucial ingredient in epidemiological studies, but also for other fields such as transportation planning.

II. STATISTICS OF INTERURBAN FLOWS

We will first study standard statistical properties of interurban flows obtained from migration data collected by Baidu Qianxi (see Material and Methods). This dataset enables us to monitor the traffic flows between cities. For each day d (d = 1, 2, . . . , T ), we can extract the number of individuals N (i, j, d) going from city i to city j with any travel mode. The migration data can thus be seen as a directed, weighted network of flows between the set of n = 296 cities of China whose populations are also known (see Material and Methods). We collected the data for the Spring Festival of 2020 (from Jan. 1st to Feb. 12th, 2020) and for assessing the impact of travel bans, we also collected the data for the Spring Festival of 2019 (which according to the Chinese lunar calendar takes place from Jan. 12th to Feb. 23rd, 2019).

Large heterogeneity of flows

We first consider the distribution of all flows of individuals N (i, j, d) for all cities i and j and all days d and which is shown in Fig. 1 (a). The maximum flow is of order 10 5 and the average of order 10 3 indicating a broad distribution. A power law fit is consistent with this picture with an exponent α ≈ 2.3 (Fig. 1 (a)). This heterogeneity is confirmed in Fig. 1 (b) which shows both the average value µ d and the standard deviation σ d computed over all inter-city flows (for each day d). We see that for most days the relative dispersion σ d /µ d is of order 5 -10. This heterogeneity is probably due to the large diversity of cities that can serve as origins or destinations of flows (see below for further analysis). An important feature that we can observe on Fig. 1 (b) is the sharp drop of the standard deviation after Jan. 25th, the Lunar New Year (LNY), which we will see below is mainly due to the travel ban (see also Fig. S1,S3 in SI for a detailed discussion).

Temporal versus spatial fluctuations

In order to understand the nature of the different fluctuations affecting the flows N (i, j, d), we compute the relative standard deviation ∆ ij = σij µij , where µ ij and σ ij are the average and standard deviation computed over µ d , average over all flows, for a given day d. We show on Fig. 2 (a) the spatial dispersion ∆ d versus time and in Fig. 2 (b) the distribution of ∆ ij . We observe that the spatial dispersion is of order 8.3, while the temporal dispersion is less (mainly concentrated around 1). The main reason for heterogeneity thus lies in the flow fluctuations between different origins and destinations, while temporal fluctuations are smaller but not negligible. These two sources of heterogeneity clearly represent a challenge for modeling these flows, especially with very simplified models.

Our results indicate that the first modeling step would be to describe the spatial heterogeneity of flows and then consider temporal variations.

The next natural quantities that can be computed over this network are the incoming and outgoing flows defined by

N in (i, d) = n j=1 N (j, i, d) N out (i, d) = n j=1 N (i, j, d) (1)
respectively. We measure in the same way as above various measures of fluctuations, either averaged over cities or over time, leading to the quantities

∆ in (out) d , ∆ in (out) i .
As these quantities are sums of random variables, we expect smaller relative dispersions than for N (i, j, d) which is indeed what we observe (see Fig. 2 (c) and (d), with typical values of relative dispersion of order 1 (see Figs. S4,S5 for additional details). In order to get first insights about the influence of travel bans, we compare the incoming flows and outgoing flows versus city population in 2019 and 2020 with days N before in, out before and N after out,in after LNY. We first observe that (see Fig. S2 in Supplementary Information (SI)) basically the number of outgoing individuals before LNY corresponds approximately to the number of incoming individuals after LNY with N before in (out) ≈ N after out (in) (and vice-versa). These relations thus correspond roughly to the conservation of the number of individuals traveling during the Chinese Spring Festival.

Structure of incoming and outgoing flows

The value of incoming or outgoing flows gives information about the volume of migrations, but not about the number of important origins or destinations. In order to characterize the dispersion over different cities, we denote by O(i, d) and D(i, d), the sets of origin of flows incoming in city i and destinations of flows from city i (for the day d), respectively. We then use Gini indices [START_REF] Philip M Dixon | Bootstrapping the gini coefficient of inequality[END_REF] that capture the dispersion of incoming and outgoing flows and are given by

G in (i, d) = 1 2O 2 N in (i, d) p,q∈O(i,d) |N (p, i, d) -N (q, i, d)| (2)
G out (i, d) = 1 2D 2 N out (i, d) p,q∈D(i,d) |N (i, p, d) -N (i, q, d)| (3)
where O and D represent the number of elements of the sets O(i, d) and

D(i, d). The quantity N in (i, d) = Nin(i,d) |O(i,d)|
is the average incoming flows and

N out (i, d) = Nout(i,d) |D(i,d)|
the average outgoing flows. Intuitively, if all traffic flows to city i are from one single origin city on day d, the Gini index G in (i, d) will be 1, while if traffic flows to city i are all equal, the Gini index G in (i, d) will be 0 (and similarly for G out (i, d)).

We plot these Gini indices computed for each city versus the traffic flows to or from this city. These figures 3 (a,b) show that on average the larger the traffic flows are, the more dispersed they are over a larger number of origins or destinations. In terms of epidemic control, it is clear that cities with a large flow N in and a small Gini index G in are the most critical, in the sense that many people from many different cities are converging to the same place. Equally, cities with a large N out and a small G out should be particularly monitored, since they can act as hubs in spreading the disease over the inter-city network. As shown in Fig. 3 (c,d), we show the top 5 critical cities, including Beijing, Shanghai, Chongqing and Guangzhou for both the incoming and outgoing flows, Shenzhen for the incoming flows, and Dongguan for the outgoing flows.

III. STATISTICAL STRUCTURE OF THE NATIONAL POPULATION

An important effect of incoming and outgoing flows is that they change the population structure. Some cities will receive a large number of individuals while for others we expect a decrease of their population. Migration thus affects the statistical structure of the national population and in this section we will characterize this effect.

Temporal evolution of population structure

In order to characterize the disparity of the population distribution and how it varies during seasonal migrations, we consider the population of city i at time d given by

P (i, d) = P 0 (i) + d d N in (i, d ) - d d N out (i, d ), (4)
where P 0 (i) represents the population of city i without incoming and outgoing flows. The Gini index for the city population of the whole country at day d is then given by

G(d) = 1 2n 2 P (d) n i,j=1 |P (i, d) -P (j, d)|, (5)
where

P (d) = 1 n n i=1 P (i, d)
is the average population of all cities at day d. Intuitively, if all people gather in one city, G will be 1, while if people spread evenly across all cities, G will be 0. For comparison, we also define the Gini index at rest as

G rest = 1 2n 2 P 0 n i,j=1 |P 0 (i) -P 0 (j)|. (6)
This quantity captures the degree of population concentration without any traffic flows, where P 0 = 1 n n i=1 P 0 (i) is the average population of all cities without any traffic flows. We show in Fig. 4 the variation of the Gini coefficient when we take into account migration flows. We plot both the results for 2019 and 2020. In both cases we see an important increase of the Gini coefficient in a short time (about a week): when the LNY is approaching, people go back from workplaces to hometowns for reunion with families. A smaller set of cities concentrates these meetings with the number of important cities reaching its minimum and the Gini index reaching its peak on the LNY. Based on the Gini index, we can estimate the number of 'important' cities where the concentration takes place through [n(1 -G)], where [•] denotes the integer part (see SI for details where we show in Fig. S6 this number versus time and we indeed observe an important drop when approaching the LNY). After the LNY (Jan. 25th), individuals are going back home and the Gini coefficient relaxes back to its original value, but much slower. We observe that in 2020, the increase of the Gini index is larger and, due to travel bans, the decrease even slower than normal. The reason may be that after the outbreak of COVID-19, almost all regions have deferred the time of resuming works and classes after the Spring Festival holiday. For example, Shanghai proposed that companies not crucial to the nation should not resume works before Feb. 10th and that schools should provide online classes. At this point, the population structure at the national level is far from being back to normal. These different results show that these seasonal movements induce a strong concentration of individuals in a relative small set of cities, and that travel bans tend to keep this situation of high concentration.

Return to 'equilibrium': pendular ratio

We observe in Fig. 4 that after the LNY there is a decrease of the Gini index indicating a return to normal state characterized by a lower concentration of individuals. In order to characterize quantitatively this return to the original state (before holidays), we measure the gap between individuals going out from a city before the LNY and coming back after it. This gap defines a 'pendular ratio' given by

R(i, d f ) = d<d d+d f N in (i, d) d-d f d< d N out (i, d) , (7)
where d f is a range of days around the LNY d. If this ratio is much larger than 1, it means that for this city there is a large incoming flow while for the opposite situation R(i, d f ) 1, a large number of individuals are going out (compared to the incoming flows). At large times d f , we expect that R 1 since most of the in- dividuals have come back. We divide cities into three categories according to the value of R(i, 1): If the value is larger than 1.5, we classify city i as a 'receiver' city. If the value is less than 0.5, we classify city i as an 'emitter' city. Finally, if the value is between 0.5 and 1.5, we classify city i as a 'transit' city. We represent on Fig. 5 the cities of different types on the map of China. We observe that both receiver and transit cities are homogeneously distributed in China. In constrat emitters cities are in general located in developed regions, e.g., Beijing, Shanghai, Guangzhou, and so on, as shown in figures 5 (a) and (b). It is interesting to note that cities of the Hubei province (within the dashed circle in the figure) are emitters cities in 2020, essentially due to travel restrictions that prevented individuals to come back to Wuhan. This is an important difference compared to the year of 2019 that appears here in the spatial structure of emitters and receivers. We show in Fig. 6 (a,b) the pendular ratio for 2019 and 2020 for all cities and we highlight 5 cities: Wuhan, Beijing, Tianjin, Chongqing and Shanghai, corresponding to the origin place of COVID-19 and four province-level municipalities. We note here that the curve corresponding to Wuhan is at the bottom of all cities in Fig. 6 (b), reflecting the success of sealing off Wuhan from all outside contact to stop the spread of the disease since Jan. 23rd.

In Fig. 6 (c,d) we show this pendular ratio for 2019 and 2020 for the different types of cities (we average over cities in a given category, emitter, receiver or transit). We observe that the standard deviation is small for the three groups adding credit to their definition. In addition, compared to 2019, the values of R(i, 1) corresponding to 2020 are much smaller. We observe that in 2019, the pendular ratio of all the three types of cities returns to 1, meaning that the majority of individuals who went away for the holidays came back. The situation for 2020 is very different with a pendular ratio for all types of cities that converges to a value less than 1 (even less than 0.5), indicating that the majority of people who went away for the holidays did not come back yet. This result remains consistent with the conclusion of Gini index (Fig. 3) about a larger concentration in cities and the effect of travel bans.

Finally, we note here that we additionally implemented our whole analysis at the province level (see SI) and the results obtained are similar are those obtained at the city level.

IV. CONCLUSION

Our findings thus concern four different aspects. First, the traffic flows between cities are very heterogeneous not only spatially but also from a temporal perspective. Such a large heterogeneity could be induced by both the Spring Festival and the travel ban. Similar results apply also to an aggregated level, i.e. the incoming and outgoing flows for cities also display important heterogeneities. This aspect is crucial for understanding and modeling epidemic spreading for which we know the importance of heterogeneity for epidemic spreading on networks [START_REF] Pastor | Epidemic spreading in scale-free networks[END_REF][START_REF] Barthélemy | Velocity and hierarchical spread of epidemic outbreaks in scale-free networks[END_REF] and more generally for most processes [START_REF] Barrat | Dynamical processes on complex networks[END_REF]. We also quantify the dispersion of origins/destinations of the incoming/outgoing flows showing that for larger flows we have a larger variety of origins and destinations. We also show that during these seasonal migrations of the Spring Festival, the national structure of population changes quickly with a larger concentration in a small set of cities. This concentration decays normally in time after the festivities but travel bans slow down this return to the initial state. It is natural to try to stop the geographical spread of the disease by stopping interurban movements, but on the other hand, large concentration in cities can favor the spread at the city level and increase the number of infected cases. This concentration can be compensated by a more important control at the individual contact level which is what was done in cities such as Wuhan. These results are in line with epidemic modeling results [START_REF] Chinazzi | Preliminary assessment of the international spreading risk associated with the 2019 novel coronavirus (2019-ncov) outbreak in wuhan city[END_REF], where it was shown that travel quarantine is effective only when combined with a large reduction of intracommunity transmission. Our results thus highlight the importance of mobility studies for modeling a variety of processes and in particular for understanding and modeling the spread of epidemics. Effective mitigating strategies need to take into account the change of population structure that we exhibited here.

V. MATERIAL AND METHODS

Data

We obtained the migration data from Baidu Qianxi (http://qianxi.baidu.com), based on Baidu Location Based Services and Baidu Tianyan, for all transportation modes. It provides the following two datasets: migration index reflecting the size of the population moving into or out from a city/province, and migration ratio capturing the proportion of each origins and destination. We collected the data during Chinese Spring Festival period of 2020 (from Jan. 1st to Feb. 12th, 2020). For parallel comparison, the migration index during the same period of 2019 (re-scaled according to Chinese lunar calendar, from Jan. 12th to Feb. 23rd, 2019) is also used.

In addition to the migration data, we collected the demographic from China Statistical Yearbook (http://www.statsdatabank.com), an annual statistical publication, which reflects comprehensively economic and social development of China. It covers key statistical data in recent years at both the city level and the province level. We collected the data of population of 31 province-level regions and 296 city-level regions from China Statistical Yearbook 2019, the latest edition provided.

deviations corresponding to incoming and outgoing flows with smaller relative dispersions are shown in figures S10 (c) and (d).

We compare the incoming flows and outgoing flows versus city population in 2019 and 2020 at province-level with days before and after LNY highlighted by different colors in Fig. S11. Compared to 2019 (figures S11 (a) and (b)), the differences between days before and after LNY are much larger in 2020 (figures S11 (c) and (d)).

Statistical structure of the national population

The trends of the population Gini index for 2019 and 2020 are shown in Fig. S12. The Gini index reaches its maximum around LNY and returns to normal state gradually. Compared to 2019, the Gini index corresponding to 2020 has a higher peak and decreases with a slower speed.

We apply the criteria of three categories of provinces similar to that for city: If the value of R(i, 1) is larger than 1.2, we classify province i as a 'receiver' province. If the value is less than 0.8, we classify province i as an 'emitter' province. Finally, if the value is between 0.8 and 1.2, we classify province i as a 'transit' province. We show the location of three categories of cites on map of China for 2019 and 2020 in Fig. S13 and observe that receiver provinces are the majority in 2019 while emitter provinces are the majority in 2020. This results seem to make sense since most people defer the return time due to the travel ban, so that most provinces are 'emitters' in 2020.

We show in figures S14 (a) and (b) the pendular ratio for 2019 and 2020 for all provinces and highlight 5 provinces: Hubei, Beijing, Tianjin, Chongqing and Shanghai, corresponding to the origin province of COVID-19 and four province-level municipalities. We note that the curve corresponding to Hubei is in the bottom from all provinces in Fig. S14 (b), indicating that except Wuhan, the origin city of COVID-19, people also avoid going to cities of Hubei. We observe that, in 2019, the pendular ratios of all the three types of cities return to 1, meaning that the majority of individuals who went away for the holidays came back, as shown in Fig. S14 (c). The situation for 2020 is very different with a pendular ratio for all types of cities that converges to a value less than 1, indicating that the majority of people who went away for the holidays did not back yet, as shown in Fig. S14 (d).

To sum up, we show that the traffic flows between provinces are very heterogeneous, and display both large temporal and spatial fluctuations, and so on. These results for province-level are in good agreement with for city-level, indicating that our methods are applicable to both scales. Despite the detailed characters of traffic flows revealed by results for city-level and province-level, a global view of a higher level is also necessary. The statistical properties of the interurban mobility help us to understand the effect of travel restrictions, their impact on and the control of epidemic spread.

Figure 1 :

 1 Figure 1: (a) Distribution of all traffic flows N (i, j, d) in loglog. The line is a power law fit of the form N -α with exponent α = 2.27 with fitting method described in [26]. (b) Average and standard deviation of the flows N (i, j, d) averaged over traffic flows versus the date d (from 1st January to 12th February).

Figure 2 :

 2 Figure 2: (a) Relative standard deviation of the flows N averaged over traffic flows and represented here versus time. (b) Distribution of the relative standard deviation of N averaged over time. (c) Relative standard deviation of Nin(i, d) and Nout(i, d) averaged over cities and shown here versus time. (d) Distribution of the relative standard deviation of Nin and Nout averaged over time

Figure 3 :

 3 Figure 3: (a) Gin versus Nin for all cities and days. (b) Gout versus Nout for all cities and days (shown in loglog). The thick line indicates the average of Gin ( Gout ) versus of Nin (Nout) and the shaded area represents the standard deviation of the average. Top 5 critical cities for (c) incoming flows and (d) outgoing flows.

Figure 4 :

 4 Figure 4: Temporal variations around the Spring Festival holidays of the population Gini index for 2019 and 2020. The horizontal dotted line represents the value 'at rest'. The vertical line indicates the day of the LNY.

Figure 5 :

 5 Figure 5: (a) Emitter, receiver and transit cities according to the value of R(i, 1) for 2019 with the number of three categories of cites at lower left corner. (b) Emitter, receiver and transit cities according to the value of R(i, 1) for 2020 with the number of three categories of cites at lower left corner.

Figure 6 :

 6 Figure 6: Comparison of pendular ratio between 2019 and 2020. (a) Pendular ratio for all cities versus days d f from LNY in 2019. We highlight five important cities. (b) Pendular ratio for all cities versus days d f from LNY in 2020 with highlight of 5 cities. (c) Average values of the pendular ratio over cities according to the classification (receiver, emitter or transit cities) versus days from LNY in 2019. The colored areas correspond to one standard deviation. (d) Mean value of the pendular ratio over cities according to the classification versus days from LNY in 2020 with the shaded area representing the corresponding standard deviation.

Figure S1 :

 S1 Figure S1: Average and standard deviation of N over traffic flows versus time for 2019 with the corresponding migration ratio for 2020.

Figure S2 :

 S2 Figure S2: Incoming flows (a) and outgoing flows (b) for all cities and days versus city population before and after the 2019 LNY in loglog. Incoming flows (c) and outgoing flows (d) for all cities and days versus city population before and after the 2020 LNY in loglog.

Figure S3 :

 S3 Figure S3: Observation of incoming and outgoing flows for 2020. (a) Distribution of all incoming flows in loglog with parameters of power law fitting on the top middle. (b) Distribution of all outgoing flows in loglog with parameters of power law fitting on the top middle. (c) Average and standard deviation of Nin and Nout over cities versus time.

Figure S4 :

 S4 Figure S4: Observation of incoming and outgoing flows for 2019. (a) Distribution of all incoming flows in loglog with parameters of power law fitting on the top middle. (b) Distribution of all outgoing flows in loglog with parameters of power law fitting on the top middle. (c) Average and standard deviation of Nin and Nout over cities versus time. (d) Relative standard deviation of Nin and Nout over cities versus time. (e) Distribution of the relative standard deviation of Nin and Nout over time.

Figure

  Figure S5: (a) Gin versus Nin for some cities and all days. (b) Gout versus Nout for some cities and all days.

Figure

  Figure S6: (a) Distribution of populations of all cities in loglog. (b) The number of important cities versus time corresponding to 2019 and 2020. We match the time scale for 2019 and 2020 according to LNY (for the sake of clarity, we show on the x-axis the dates for 2020 only). The vertical line highlights the LNY.

Figure

  Figure S7: (a) Emitter, receiver and transit cities according to the value of R(i, 1) for 2019 with the number of three categories of provinces at the lower left corner. (b) Emitter, receiver and transit cities according to the value of R(i, 1) for 2020 with the number of three categories of provinces at the lower left corner.

Figure

  Figure S8: (a) Mean value of the pendular ratios over cities according to the classification (receiver, emitter or transit cities) with different criteria versus days from LNY in 2019. The colored areas correspond to one standard deviation. (b) Mean value of the pendular ratios over cities according to the classification with different criteria versus days from LNY in 2020 with shaded areas representing standard deviation.

Figure

  Figure S9: (a) Distribution of all traffic flows N (i, j, d) in loglog. The line is a power law fit of the form N -α with exponent α = 2.27. (b) Average and standard deviation of N over traffic flows versus time.

Figure

  Figure S10: (a) Relative standard deviation of N over traffic flows versus time. (b) Distribution of the relative standard deviation of N over time. (c) Relative standard deviation of Nin and Nout over provinces versus time. (d) Distribution of the relative standard deviation of Nin and Nout over time.

Figure S11 :

 S11 Figure S11: Incoming flows (a) and outgoing flows (b) for all provinces and days versus province population before and after the 2019 LNY in loglog. Incoming flows (c) and outgoing flows (d) for all provinces and days versus province population before and after the 2020 LNY in loglog.

Figure S12 :

 S12 Figure S12: Temporal variations during the Spring Festical of the population Gini index for 2019 and 2020. The dotted line represents the value 'at rest'.

Figure

  Figure S13: (a) Emitter, receiver and transit provinces according to the value of R(i, 1) for 2019 with the number of three categories of provinces at lower left corner. (b) Emitter, receiver and transit provinces according to the value of R(i, 1) for 2020 with the number of three categories of provinces at lower left corner.

Figure S14 :

 S14 Figure S14: Comparison of pendular ratio between 2019 and 2020. (a) Pendular ratio for all provinces versus days d f from LNY in 2019. We highlighted five important provinces. (b) Pendular ratio for all provinces versus days d f from LNY in 2020 with highlight of 5 provinces. (c) Mean values of the pendular ratio over provinces according to the classification (receiver, emitter or transit provinces) versus days from LNY in 2019. The colored areas correspond to one standard deviation. (d) Mean value of the pendular ratio over provinces according to the classification versus days from LNY in 2020 with colorbar representing standard deviation.

Appendix A: Supplementary Materials

Data for 2019

In order to evaluate the heterogeneity of flows of 2019 with comparison to that of 2020, we use the migration index during the same period of 2019 (re-scaled according to Chinese lunar calendar, from Jan. 12th to Feb. 23rd, 2019), and we would compute the distribution of all flows N (i, j, d), for all cities i and j and all days d, though N (i, j, d) = N out (i, d) × p(i, j, d). The Chinese Lunar New Year of 2019 is Feb. 5th. Here, N out (i, d) is migration index reflecting the size of the population moving into or out from a city/province, and p(i, j, d) is migration ratio capturing the proportion of each origins and destination. However, the migration ratio is unavailable for 2019. We apply the data of p(i, j, d) for 2020 to the computation of N (i, j, d) for 2019, with results shown in Fig. S1. This result exhibits large heterogeneity of flows and displays a localized drop around LNY.

Statistics of N in and Nout

Versus population

We observe a power law relationship between the incoming/outgoing flows and city population in Fig. S2 which indicates that the larger a city and the more flows it carries. Compared to 2019, the differences between the scatter points for incoming flows corresponding to days before and after LNY are much larger in 2020. This result emphasizes again that travel ban causes indeed the sharp drop of standard deviation in Fig. 1 (b) of the main text rather than the low travel intention during the Spring Festival.

The power law fits that we obtain imply that N in ∼ P γin 0 where γ in ≈ 0.93 before LNY and γ in ≈ 0.88 after LNY. Similarly for outgoing flows we obtain γ out ≈ 0.85 before and γ out ≈ 0.93 after LNY. We can interpret these results as a consequence of the conservation of the number of individuals traveling before and after LNY.

Distribution of Nin,out for 2020

We show the distribution, average and standard deviation of incoming/outgoing flows in Fig. S3. We observe that these distributions are relatively broad, in particular outgoing flows (Fig. S3 (a) and (b)). We show the standard deviation of incoming flows and outgoing flows over cities for each day, and the corresponding average (the same for the incoming and outgoing flows) in Fig. S3 (c). Note that the standard deviation of N in is smaller than that of N out before Jan. 25th, the 2020 LNY, while the situation reverses after Jan. 25th. A reasonable explanation for this is that people go to a relative large number of hometowns from a relative small number of workplaces before the Spring Festival; Due to the travel ban, people do not come back after the Spring Festival.

Distribution of Nin,out for 2019

We show the same quantities as above but for the year 2019. Here also, we observe broad distributions both for all incoming flows and all outgoing flows are shown in Fig. S4.

Fluctuations are larger around LNY where the total flow of individuals is larger, allowing for more heterogeneity. Before LNY individuals move from a large variety of cities to a relatively small number of hometowns explaining the large fluctuations of N out . After LNY, individuals are returning from a small number of hometowns to a large variety of cities, inducing large fluctuations of N in . The corresponding dispersion and relative dispersions ∆ out d ) and ∆ in d ) are shown in figures S4 (c) and (d). This also results in the heterogeneous distribution for the relative standard deviation of N in and N out averaged over cities for 2019 in Fig. S4 (e).

Gini indices for Nin,out

We compute Gini indices for cities. Instead of showing results for all cities, we plot Gini indices versus the traffic flows to or from cities in set {i ∈ V| min

5} and Wuhan in Fig. S5. In this case from these cities, many people go out to or come in from many different cities. These cities are critical and include Shanghai, Beijing, and so on. Due to travel bans, Wuhan exhibits specific features of scatter points with clear separation before and after LNY (see the scatter points in blue at the upper left corner corresponding to Wuhan in Fig. S5).

We also observe here in both cases a decreasing behavior on average. This is more salient for N out where the trend is clearly visible. This indicates that for larger outgoing flows, the Gini is smaller with no clearly dominant flow.

Statistical structure of the national population

We first show the population distribution in Fig. S6 (a). We observe a broad distribution and a power law fit gives the exponent α ≈ 5. We also show the number of important cities quantified by the integer part of n[1 -G(d)], in Fig. S6 (b).

Pendular ratio

In order to test the dependence of the pendular ratios on the criteria for defining classes, we change the criteria as follows: here if the value of R(i, 1) is larger than 1.2, we classify city i as a 'receiver' city. If the value is less than 0.8, we classify city i as an 'emitter' city. Finally, if the value is between 0.8 and 1.2, we classify city i as a 'transit' city. We show the location of three categories of cites on the map of China for 2019 and 2020 in Fig. S7.

Compared to the criteria in the main text, the number of transit cities decreases, while the number of emitter and receiver cities increases. However, as shown in Fig. S8, the patterns of the average value of pendular ratio corresponding to three categories of cites remain unchanged.

Statistics at the inter-province level

Statistics of flows

In what follows, we show some corresponding results based on province-level data instead of city-level data. The distribution of traffic flows between provinces exhibits large heterogeneity with the exponent of power law fit α around 2.24, as shown in Fig. S9 (a). A sharp drop of the standard deviation after Jan. 25th is observed in Fig. S9 (b).

We show the relative standard deviation of N over flows versus time with an order around 2.96 in Fig. S10 (a) and the distribution of the relative standard deviation of N over time concentrating around 0.64 in Fig. S10 (b). Large heterogeneity of traffic flows between provinces confirms the difficulty of modeling these flows. The relative standard


Sylvia G Lehmann
email: sylvia.lehmann@univ-grenoble-alpes.fr

Benjam Gilbert
email: bgilbert@lbl.gov

Thie Gg Maffeis
email: t.g.g.maffeis@swansea.ac.uk

Alexei Grichine
email: alexei.grichine@univ-grenoble-alpes.fr

Isabelle Pignot-Paintrand
email: isabelle.paintrand@grenoble-inp.fr

Simon Clavaguera
email: simon.clavaguera@cea.fr

Walid Rachidi
email: walid.rachidi@cea.fr

Michel Seve
email: michel.seve@univ-grenoble-alpes.fr

Laurent Charlet
email: charlet38@gmail.com

In Vitro Dermal Safety Assessment of Silver Nanowires after Acute Exposure: Tissue vs. Cell Models

Keywords: silver nanowires, cytotoxicity, skin irritation in vitro, primary keratinocytes, 3D reconstructed epidermis model

Silver nanowires (AgNW) are attractive materials that are anticipated to be incorporated into numerous consumer products such as textiles, touchscreen display, and medical devices that could be in direct contact with skin. There are very few studies on the cellular toxicity of AgNW and no studies that have specifically evaluated the potential toxicity from dermal exposure. To address this question, we investigated the dermal toxicity after acute exposure of polymer-coated AgNW with two sizes using two models, human primary keratinocytes and human reconstructed epidermis. In keratinocytes, AgNW are rapidly and massively internalized inside cells leading to dose-dependent cytotoxicity that was not due to Ag + release. Analysing our data with different dose metrics, we propose that the number of NW is the most appropriate dose-metric for studies of AgNW toxicity. In reconstructed epidermis, the results of a standard in vitro skin irritation assay classified AgNW as non-irritant to skin and we found no evidence of penetration into the deeper layer of the epidermis. The findings show that healthy and intact epidermis provides an effective barrier for AgNW, although the study does not address potential transport through follicles or injured skin. The combined cell and tissue model approach used here is likely to provide an important methodology for assessing the risks for skin exposure to AgNW from consumer products.

Introduction

In the past 10 years, silver nanowires (AgNW) have emerged as an attractive material for a range of applications including transparent, flexible, conductive electrodes and thin films. With competitive properties to indium tin oxide (ITO), such as high optical transparency and low sheet resistance, but at lower cost and with resistance to flexing, they are predicted to replace ITO in solar cells, touch screen, liquid and flexible displays [START_REF] Liu | Silver nanowire-based transparent, flexible, and conductive thin film[END_REF][START_REF] Hu | Scalable coating and properties of transparent, flexible, silver nanowire electrodes[END_REF][START_REF] Langley | Flexible transparent conductive materials based on silver nanowire networks: A review[END_REF].

Although many studies have investigated the toxicity of silver nanoparticles (AgNP), finding a significant role for the release of Ag + ions and oxidative stress induction [START_REF] Jennifer | Nanoparticle Technology as a Double-Edged Sword: Cytotoxic, Genotoxic and Epigenetic Effects on Living Cells[END_REF][START_REF] Dubey | Perturbation of cellular mechanistic system by silver nanoparticle toxicity: Cytotoxic, genotoxic and epigenetic potentials[END_REF], few studies have been dedicated to AgNW. Key studies using models for respiratory exposure have demonstrated that AgNW confirm at least partially to the fibre pathogenicity paradigm, with a length threshold for pulmonary inflammation and toxicity in rodents and frustrated phagocytosis by macrophages in vivo and in vitro [START_REF] Schinwald | Use of silver nanowires to determine thresholds for fibre length-dependent pulmonary inflammation and inhibition of macrophage migration in vitro[END_REF][START_REF] Schinwald | Use of back-scatter electron signals to visualise cell/nanowires interactions in vitro and in vivo; frustrated phagocytosis of long fibres in macrophages and compartmentalisation in mesothelial cells in vivo[END_REF][START_REF] Schinwald | The threshold length for fiber-induced acute pleural inflammation: Shedding light on the early events in asbestos-induced mesothelioma[END_REF][START_REF] Silva | Short versus long silver nanowires: A comparison of in vivo pulmonary effects post instillation[END_REF][START_REF] Chung | Inactivation, Clearance, and Functional Effects of Lung-Instilled Short and Long Silver Nanowires in Rats[END_REF]. Additional in vitro studies with alveolar epithelial cells, non-immune cell models for respiratory exposure, have observed variable measures of toxicity. In A549 cells, exposure to AgNW induces some cytotoxicity in one study [START_REF] Stoehr | Shape matters: Effects of silver nanospheres and wires on human alveolar epithelial cells[END_REF], but not in another one at low doses [START_REF] Singh | Silver nanowires as prospective carriers for drug delivery in cancer treatment: An in vitro biocompatibility study on lung adenocarcinoma cells and fibroblasts[END_REF]. In other cell lines, such as TT1 and ATII, no toxicity was observed at the concentration tested although an increased inflammatory response was observed [START_REF] Chen | Sulfidation of silver nanowires inside human alveolar epithelial cells: A potential detoxification mechanism[END_REF][START_REF] Sweeney | Silver nanowire interactions with primary human alveolar type-II epithelial cell secretions: Contrasting bioreactivity with human alveolar type-I and type-II epithelial cells[END_REF].

Although respiratory exposure is certainly relevant to occupational health, dermal exposure may be a more important risk for consumers due to the development of touchscreen displays [START_REF] Cann | High performance transparent multi-touch sensors based on silver nanowires[END_REF][START_REF] Sannicolo | Metallic Nanowire-Based Transparent Electrodes for Next Generation Flexible Devices: A Review[END_REF], wearable textiles [START_REF] Doganay | Silver nanowire decorated heatable textiles[END_REF][START_REF] Hsu | Personal thermal management by metallic nanowire-coated textile[END_REF][START_REF] Nateghi | Silver nanowire-functionalized cotton fabric[END_REF], and medical devices incorporating AgNW [START_REF] Myers | Wearable silver nanowire dry electrodes for electrophysiological sensing[END_REF][START_REF] Jiang | Highly Efficient, and Fully Flexible Circularly Polarized Antenna Enabled by Silver Nanowires for Wireless Body-Area Networks[END_REF]. Anticipating the development of AgNW-enabled touchscreens, Verma et al. exposed a range of cell lines (epithelial, endothelial, gastric, and phagocytic) to 3-6 µm long AgNWs, in suspension or embedded in a transparent electrode, finding for the suspension AgNW a low level of cytotoxicity at a relatively low maximum dose (5 µg/mL) varying depending on cell type and increasing with AgNW doses and time of incubation [START_REF] Verma | Autophagy induction by silver nanowires: A new aspect in the biocompatibility assessment of nanocomposite thin films[END_REF]. These studies point to risks for AgNW exposure to numerous cell lines but have not provided a clear understanding of the conditions leading to adverse outcomes and are just an initial step for evaluating the potential hazard for the consumer.

Here we describe an investigation into the dermal route of exposure to AgNW. We performed paired AgNW uptake and cytotoxicity studies in human primary keratinocytes (HPK) cells and in reconstructed human epidermis (RHE) model (epiCS). Indeed, HPK are the major cell component of the epidermis and RHE is an alternative model that presents the three-dimensional architecture of skin tissue, including the outer protective barrier. RHE models are used for dermal exposure to cosmetics and environmental irritants and pollutants with approved OECD guidelines like the in vitro skin irritation assay [START_REF]Test No. 439: In Vitro Skin Irritation: Reconstructed Human Epidermis Test Method[END_REF]. Such an approach has been used to determine the extent of penetration and the toxicity of several nanomaterials [START_REF] Vinardell | In Vitro Comparative Skin Irritation Induced by Nano and Non-Nano Zinc Oxide[END_REF][START_REF] Kishore | Assessment of the dermal and ocular irritation potential of multi-walled carbon nanotubes by using in vitro and in vivo methods[END_REF], but not yet for AgNW.

Materials and Methods

Silver Nanowires Characterization

Short and long polyvinylpyrrolidone (PVP) coated AgNW suspensions (S-AgNW and L-AgNW) were purchased from nanoComposix (San Diego, CA, USA) and were stored in lightproof bottles in an anaerobic chamber. Both solutions of AgNW were imaged by SEM and the length and diameter observed were in agreement with manufacturer specifications and with previous observations [START_REF] Scanlan | Silver nanowire exposure results in internalization and toxicity to Daphnia magna[END_REF] (Table S1).

AgNW were supplied into water at 1 mg/mL suspension. They were diluted into complete keratinocytes culture medium, then agitated gently for 1 min and they were well dispersed without aggregation, as observed throughout this study by light microscopy (data not shown). AgNW settling rates were assessed by measuring a series of optical absorption spectra (Ocean Optics, Largo, FL, USA) of unstirred 4-mL suspensions of 1 µg/mL S-AgNW and L-AgNW in keratinocyte culture media (KSFM; Life Tech., Villebon Sur Yvette, France) at different times. AgNW present a peak of absorption at 380 nm, the decrease of absorbance at 380 nm over the time reflect then AgNW settling. Measurements were performed in closed 1-cm path-length quartz cuvettes.

The calculation of AgNW concentration expressed as the number of AgNW/mL was made using the following formula:

N = C/(π × r 2 × l × d)
where N = number of wires/mL; C = concentration (g/mL); π = Pi; r = radius; l = length; d = density. The calculated volume and surface area/single NW is indicated in Table S1.

Ionic Silver Release into Media

Ag + release into complete keratinocytes culture medium was quantified by adding AgNWs to culture medium and measuring the dissolved silver concentration at different intervals with ICP-MS (Perkin-Elmer Elan DRC II) as previously described [START_REF] Scanlan | Silver nanowire exposure results in internalization and toxicity to Daphnia magna[END_REF].

Cell Culture

Normal human keratinocytes were isolated from skin tissue obtained after breast plastic surgery, from healthy Caucasian donors (32 to 34 year old) with their informed consent (Department of "Chirurgie Plastique et Maxillo-faciale", CHU Grenoble, La Tronche, France). Keratinocytes were isolated from the skin explant as described previously [START_REF] Mouret | Cyclobutane pyrimidine dimers are predominant DNA lesions in whole human skin exposed to UVA radiation[END_REF] and cells from at least three different donors were used for each experiment. Human primary keratinocytes were then cultured in keratinocyte serum free medium (KSFM; Life Tech.) supplemented with 25 µg/mL Bovine Pituitary extract (BPE), 0.9 ng/mL recombinant human epithelial growth factor (EGF) and primocin (Life Tech., Villebon Sur Yvette, France) at 37 • C in a humidified atmosphere containing 5% CO 2 . For all experiments cells were used at passages 2 or 3 only.

Three dimensional human reconstructed epidermis epiCS (CellSystems, GmbH, Troisdorf, Germany) were maintained according to the suppliers' instructions in 6-well plates at 37 • C, 5% CO 2, and 95% relative humidity with epiCS medium.

Cell Viability Assays

Cytotoxicity was assessed using the 3-[4,5-dimethylthiazol-2-yl]-2,5-diphenyl tetrazolium bromide (MTT) assay (Sigma-Aldrich, St. Louis, MO, USA).

Human Primary Keratinocyte (HPK) assays were performed as follows. Cells were seeded in 96-well plates in complete KSFM medium. After 24 h, cells were treated with concentrations of 1.56, 3.13, 6.25, 12.5, 25, and 50 µg/mL of L-AgNW, S-AgNW and AgNO 3 as a control for silver ion toxicity. Silver nitrate (99.999% trace metal basis) was obtained from Sigma-Aldrich. After 24 h, cells were washed two times with phosphate buffered saline (PBS) and fresh KSFM medium was added supplemented with 0.5 mg/mL MTT. After 3 h of incubation at 37 • C, formazan crystals were dissolved with 100 µL Dimethyl sulfoxide (DMSO) and then transferred into new plates for absorbance reading at 560 nm using a microplate reader (Varioskan ® Flash, Thermo Scientific, Villebon-sur-Yvette, France). Cell viability was expressed relatively of absorbance of negative control (non-treated cells). Data analyses were performed with Excel 2013 for Windows. In order to ensure that there was no interference between the MTT and AgNW, the same experiment was performed without cells and the NW were pelleted by centrifugation (15 min at 800 g) before the addition of MTT. The absorbance measured was similar to the background (data not shown).

In order to determine if Ag + ions released from the AgNW had any contribution to the observed loss of viability, we performed an ion-release experiment. After 24-h incubation with cells, cell supernatants were harvested, centrifuged (15 min at 800 g), and the supernatants were transferred onto other cells for another 24 h prior to a second MTT assay that was performed as described previously.

In Vitro Skin Irritation Assay with epiCS RHE

Reconstructed Human Epidermis (RHE) assays were performed as follows. After 24 h of equilibration time, RHE were exposed topically to 50 µL of 0.6 mg/mL of S-AgNW, L-AgNW and the equivalent concentration of Ag + as AgNO 3 . A positive control with 3% Sodium dodecyl sulfate (SDS) was also performed. All dilutions were performed in sterile water and the suspensions fully covered the RHE surface. The same amount of water was added topically to the control non-treated samples. The treatment was applied as recommended in the OECD Test Guideline 439 [START_REF]Test No. 439: In Vitro Skin Irritation: Reconstructed Human Epidermis Test Method[END_REF] (in vitro skin irritation assay) for 20 min, then the RHE was washed 20 times in PBS and incubated 42 h. We also performed a longer exposition time, when no washes step was performed and incubation lasted 42 h. At the end of the 42 h, RHE were rinsed with PBS and MTT assay was performed according to the manufacturer's instructions. Absorbance measurement and analysis of the data were done as previously stated.

Electron Microscopy Analysis (SEM and TEM)

HPK were grown on glass coverslips and incubated with 1.56 µg/mL of AgNW for 24 h. Then, HPK were rinsed with Hanks Balanced Salt Solution, fixed in 4% paraformaldehyde prepared in a 0.1-M sodium cacodylate buffer at pH 7.2 (Sigma-Aldrich, St. Louis, MO, USA), dehydrated in 5 min baths of increasing concentrations of ethanol (70; 80; 85; 90; 95, and 100%) and then the samples were allowed to air dry at RT. Whole cells were observed with a Hitachi S4800 scanning electron microscope (SEM) (Centre for Nanohealth, College of engineering, Swansea, UK). Most samples were imaged without electrical charging artefacts but a few cases were coated with a Cr coating of 5 nm.

Some HPK samples treated for 72 h were fixed in 2.5% glutaraldehyde in a 0.1-M sodium cacodylate buffer at pH 7.2, post-fixed in 1% osmium tetroxide and 1.5% potassium hexacyannoferrate(II) trihydrate, dehydrated with 30 min baths of increasing concentrations of ethanol (70; 95 and 3 times in 100%) and embedded in Epon resin (Delta Microscopy). Sections were cut (250 nm) with an ultramicrotome (Leica (Wetzlar, Germany), NanoID, CEA) and observed with a FEG SEM ZEISS GeminiSEM 500 (CMTC, France) using a backscatter detector. The contrast of images was inverted with Image J (NIH).

RHE samples that were exposed to AgNW for 42 h were rinsed with PBS, fixed in 2.5% glutaraldehyde in a 0.1-M sodium cacodylate buffer at pH 7.2, post-fixed with 1% osmium tetroxide for 1 h, counterstained with 0.5% uranyl acetate at pH 4 overnight, then dehydrated in graded concentrations of ethanol and embedded in Epon resin (Delta Microscopy). Ultra-thin and semi-fine sections (60 nm and 400 nm, respectively) with an ultramicrotome (Ultracut Leica) and observed with a JEOL 1200 EX transmission electron microscope (TEM) operated at 80 kV (Grenoble Institut des Neurosciences, France). Semi-thin sections were stained with toluidine blue and visualized with a phase contrast microscope to assess RHE morphology.

Confocal Microscopy

HPK were grown on Lab-Tek™ chamber slides and incubated for 24 h with 1.56 µg/mL of AgNW. DNA was stained by the addition of Hoechst 33342 (1 µg/mL in PBS, Sigma) and images were acquired with a LSM710 confocal microscope (IAB, Grenoble, France). Hoechst fluorescence was excited with 405 nm laser and detected in 420-480 nm interval in confocal descanned mode with a pinhole diameter of 1 Airy unit. The AgNW were detected in confocal reflected light mode. For time lapse analysis, HPK were incubated immediately before imaging with 1.56 µg/mL of AgNW. Three-dimensional (3D) stack images were acquired using a 40X/1.2 water immersion objective with Zen software (Carl Zeiss, Oberkochen, Germany). The top and bottom points of the cells were defined, and 29 z-stacks with a step size of 0.5 µm per time point were acquired each 2 min for 8 h with a LSM 510 Carl Zeiss confocal microscope (IAB, Grenoble, France). Image analysis was performed with Image J (NIH).

Cytoviva Microscopy

HPK were grown on glass coverslips and exposed for 72 h to 1.56 µg/mL of AgNW. After incubation, the cells were rinsed with PBS, fixed in 4% paraformaldehyde and stained with Hoechst 33342. Dark-field micrographs and hyperspectral images were acquired with a CytoViva hyperspectral microscopy system (CEA, NanoSafety Platform, Grenoble, France). The hyperspectral data were converted to text files using ENVI software and analysed in the IgorPro software using the GG Macros routines developed for spectromicroscopy. The lineshapes of all visible-to-near-infrared (V-NIR) spectra dominated the output of the lamp and the precise alignment of the condenser, making quantitative spectroscopy very difficult. However, by subtracting a spectrum from a nearby bright Ag-free biological structure, V-NIR spectra of AgNW and AgNR could be obtained.

Results and Discussion

AgNW Characterization

The two sizes of commercial PVP-coated AgNW used in this study came from the same supplier as prior studies and as described in the methods section. They possessed similar dimensions and minor nanorods (AgNR) contamination as described previously [START_REF] Silva | Short versus long silver nanowires: A comparison of in vivo pulmonary effects post instillation[END_REF][START_REF] Scanlan | Silver nanowire exposure results in internalization and toxicity to Daphnia magna[END_REF]. The average dimensions were approximately 2 µm × 40 nm (length × diameter) for the short samples (S-AgNW) and 20 µm × 50 nm for the long sample (L-AgNW) (Table S1). AgNW diluted into complete keratinocytes culture medium were well dispersed without aggregation, as observed throughout this study by light microscopy (data not shown). We followed AgNW settling in culture medium by measuring the decrease of absorbance at 380 nm (Figure 1A). L-AgNW settled faster than S-AgNW. Although the spectrophotometric data cannot easily be converted to settling rates, it is likely that the L-AgNW settled completely within 24 h in the well-plates used for toxicity assays. The release of Ag + into cell-free media was measured for 50 µg/mL AgNW dilution by ICP-MS over the time (Figure 1B). After 28 h incubation, Ag + concentration reached 0.4-0.5 µg/mL for the short-and long-NW respectively. Ag + dissolution was measured also after 2 weeks of incubation and the concentration was 0.6 µg/mL for both types of AgNW.

Human Primary Keratinocytes

Acute Cytotoxicity

The acute cytotoxicity of AgNW and Ag + to human primary keratinocytes (HPK) was assessed using the MTT assay. Both types of AgNW induced moderate toxicity (Figure 2A) with a shallow dose response curve showing 35% and 50% viability after 24 h at the highest concentration tested (50 µg/mL) for the S-and L-AgNW, respectively. Much higher cytotoxicity was observed for Ag + with a sharper dose response, zero cell viability at 12.5 µg/mL and a calculated LD50 of 4 µg/mL. Similar cytotoxicity was observed with the neutral red uptake assay (data not shown).

In order to investigate the possible contribution of Ag + release in the cell medium to the toxicity of AgNW, we harvested the cell supernatant after 24 h of exposure, centrifuged it, and added it to a fresh cell sample to perform an MTT assay after the 24-h incubation. All supernatants from AgNW exposure showed no toxicity at all (Figure 2B), consistent with the low measured Ag + release, approximately 1% w/w (Figure 2B), and showing that the cytotoxicity observed after AgNW treatment is not due to Ag + release in medium. Similar results were reported with A549 cells after 48 h of exposure [START_REF] Stoehr | Shape matters: Effects of silver nanospheres and wires on human alveolar epithelial cells[END_REF] although Ag + release in RPMI 1640 medium supplemented with 10% FBS was higher than observed here for the keratinocytes medium that does not contain serum. Verma et al. also performed this ion release control experiment and did not observe any toxicity of the AgNW supernatant in A549 cells cultivated in Ham's F12 medium containing 10% FBS [START_REF] Verma | Autophagy induction by silver nanowires: A new aspect in the biocompatibility assessment of nanocomposite thin films[END_REF].

The supernatant from the AgNO 3 exposures exhibited lower-toxicity responses in this second experiment than in the first one (Figure 3B). For example, while only ~10% viability was obtained in the first assay at 6.25 µg/mL Ag + , the supernatant from this sample gave 100% viability. This is consistent with a portion of the added Ag + being removed from the media by sorption or uptake. Several reports [START_REF] Oberdörster | Nanotoxicology: An emerging discipline evolving from studies of ultrafine particles[END_REF][START_REF] Sager | Surface area of particle administered versus mass in determining the pulmonary toxicity of ultrafine and fine carbon black: Comparison to ultrafine titanium dioxide[END_REF] have suggested that mass concentration may not be the most appropriate dose metric for nanotoxicology studies. We then expressed our results using other metrics and in Figure 3 we compare cell viability relatively to mass of NW (Figure 3A), particle number (Figure 3B) and surface area (Figure 3C). We observed an almost linear dose-response trend across both data sets with particle number (Figure 3B). Thus our data on AgNW support the conclusions of previous studies [START_REF] Lison | Paracelsus in nanotoxicology[END_REF][START_REF] Schmid | Surface area is the biologically most effective dose metric for acute nanoparticle toxicity in the lung[END_REF] suggesting that particle number is the most relevant dose metric for fibre like nanomaterials. Also, S and L-AgNW induced similar toxicity, in agreement with prior work in which 1.5 and 8 µm AgNW induce a similar toxicity on A549 cells [START_REF] Stoehr | Shape matters: Effects of silver nanospheres and wires on human alveolar epithelial cells[END_REF].

Scanning Electron Microscopy

SEM with secondary-electron (SE) detection provided images of the surface morphology of the whole HPK fixed cells while imaging with backscattered-electron detection (BSE) revealed both AgNW outside and inside cells with high contrast [START_REF] Schinwald | Use of back-scatter electron signals to visualise cell/nanowires interactions in vitro and in vivo; frustrated phagocytosis of long fibres in macrophages and compartmentalisation in mesothelial cells in vivo[END_REF]. Using the two detectors, SEM analysis of whole HPK cells revealed efficient internalization of AgNW after 24 and 72 h of exposure. For example, HPK cells exposed to a low dose of AgNW for 24 h showed many examples of complete AgNW internalization (Figure 4A-D) as well as partially internalized AgNW (Figure S1). No AgNW-like structures were present in any control non-treated cell (Figure S2).

SEM images of cross sectioned HPK cells also observed AgNW throughout the cell body (Figure 5). These electron-dense dots or rods were not present in cells that were not exposed to AgNW (Figure S3). Some AgNW were clearly localized inside endocytic vesicles but other AgNW were apparently present inside the cytoplasm at no specific location. The latter could have entered through an endocytic pathway and reached the cytoplasm through lysosome membrane damage, or may have pierced directly the cell membrane as demonstrated for some carbon nanotubes [START_REF] Lacerda | Cell-penetrating CNTs for delivery of therapeutics[END_REF]. Internalized AgNW were never present inside the nuclear compartment, even after 72 h of exposure, in accord with a previous study of AgNW [START_REF] Sweeney | Silver nanowire interactions with primary human alveolar type-II epithelial cell secretions: Contrasting bioreactivity with human alveolar type-I and type-II epithelial cells[END_REF]. These results contrast with silver nanoparticles that were identified inside the nucleus of IMR-90 and U251 human cell lines [START_REF] Asharani | Cytotoxicity and genotoxicity of silver nanoparticles in human cells[END_REF], although the silver nanoparticles had a smaller diameter (6-20 nm) than the AgNW. Internalized S-and L-AgNW exhibited significant morphological alteration, including bending and breaking, not observed in AgNW that settled on the substrate without cell interactions (Figure 4F). The surfaces of AgNW that were internalized inside cells also showed an increased roughness (Figure 4H), similar to prior SEM images of AgNW located in the hemolymph of AgNW exposed daphnia [START_REF] Scanlan | Silver nanowire exposure results in internalization and toxicity to Daphnia magna[END_REF]. Using high angle annular dark field scanning transmission electron microscopy (HAADF-STEM) at high magnification, Chen et al. demonstrated that sulphidization reactions leading to the formation of silver sulphide (Ag 2 S) occurred close to and on the surface of PVP-coated AgNW within human alveolar epithelial (TT1) cells [START_REF] Sweeney | Silver nanowire interactions with primary human alveolar type-II epithelial cell secretions: Contrasting bioreactivity with human alveolar type-I and type-II epithelial cells[END_REF]. It is likely that the surface modifications observed here represent similar oxidative sulphidation and reprecipitation reactions.

Confocal Fluorescence and Reflectance Microscopy

Confocal microscopy of HPK cells that were fixed and stained with Hoechst following 24-h L-AgNW exposure revealed a high density of AgNW and AgNR throughout the cytoplasm (Figure 6 and Movie S1). Longer, unbroken AgNW were clearly wrapped around the cell nucleus which neither contained any detectable AgNW nor AgNR. Time-lapse confocal microscopy of live cells captured an apparent internalization event that occurred during a 2 min interval, 1 h 26 min. after the addition of L-AgNW to cells (Movie S2). This result shows that the kinetic of AgNW internalization is fast although the mechanism remains unclear.

Hyperspectral Microscopy

Visible to near-infrared (V-NIR) spectra obtained from AgNW internalized by HPK cells consistently exhibited a red shift relative to AgNW that settled on the substrate (Figure 7). Quantification and interpretation of these spectral changes is challenging, particularly because the wavelength range of analysis does not overlap with the strong surface Plasmon resonance region of the optical excitation spectra. Moreover, the V-NIR response of Ag can be altered by several processes. For example, Pratsinis et al. used hyperspectral microscopy to locate AgNP internalized by murine macrophages (RAW 264.7 cells, ATCC) [START_REF] Pratsinis | Toxicity of Silver Nanoparticles in Macrophages[END_REF]. They observed a red-shift in some internalized objects attributed to aggregation. Leclerc and Wilkinson studied toxicity of AgNP and AgNO 3 to C. reinhardtii [START_REF] Leclerc | Bioaccumulation of Nanosilver by Chlamydomonas reinhardtii-nanoparticle or the free ion?[END_REF]. They observed nanoparticulate Ag within algal cells following exposure to either AgNP or silver ion that exhibited a distinct spectral profile relative to external AgNP. This was interpreted as being due to secondary precipitation of Ag (S) following Ag + diffusion into the cell. Because the SEM data are consistent with silver sulphidisation reactions, and because no internal aggregation of NW is observed, we interpret the spectral changes to be caused by surface chemical change.

Reconstructed Human Epidermis

Skin Irritation Assay and Histology

An in vitro skin irritation assay was performed with epiCS reconstructed human epidermis (RHE) models as recommended in the OECD test guideline 439 [START_REF]Test No. 439: In Vitro Skin Irritation: Reconstructed Human Epidermis Test Method[END_REF]. Any decrease in viability was observed after 20 min of contact with AgNW followed by extensive PBS washing and 42 h recovery (data not shown). This assay was also performed with 42 h of contact and, as described in Table 1, although the positive control (3% SDS) induced a drastic decrease of RHE viability, S-AgNW and Ag + treatment generated a weak but significant decrease in viability (96.6% and 91.9% of viability respectively). According to these results, AgNW could then be classified as non-irritant to skin because the viability was higher than 50%. Moreover, RHE treated with AgNW showed a very similar histology to control non-treated samples (Figure 8A-C).

A survey of the literature identified other nanomaterial toxicants that did not induce skin irritation using the OECD test guideline 439 and RHE models. For example, MWCNT at a dose of 10 mg/insert [START_REF] Kishore | Assessment of the dermal and ocular irritation potential of multi-walled carbon nanotubes by using in vitro and in vivo methods[END_REF] and zinc oxide nanoparticles at a concentration 0.5 mg/mL [START_REF] Vinardell | In Vitro Comparative Skin Irritation Induced by Nano and Non-Nano Zinc Oxide[END_REF] caused no decrease in RHE viability, even though the LD50 for the latter was 50 µg/mL in HaCat cell line. This difference in sensitivity between 2D cell culture and 3D RHE model is consistent with the efficient barrier protection of the reconstructed epidermis. In a recent study [START_REF] Miyani | Assessment of the in vitro dermal irritation potential of cerium, silver, and titanium nanoparticles in a human skin equivalent model[END_REF], silver nanoparticles with different sizes (10-100 nm) and with different coating (citrate, PVP or silica) as well as silver nitrate similarly did not induce skin irritation. In this study, the exposure time was 1 h followed by 42 h recovery and silver concentration was 1 mg/mL.

Transmission Electron Microscopy

We looked for evidence of cellular internalization of AgNW inside the epidermis model using TEM analysis of RHE ultrafine sections. Only one section out of five analysed showed any object consistent with the AgNW, in this case a bent NW lying on the surface of the stratum corneum layer (Figure 8D,E). Supporting elemental composition analysis was not available but the measured average diameter was in concordance with AgNW characteristics. Absence of penetration into epidermis of silver nanoparticles was also observed in vivo in pig even after topical repetitive daily application for two weeks [START_REF] Samberg | Evaluation of silver nanoparticle toxicity in skin in vivo and keratinocytes in vitro[END_REF], although, focal inflammation and oedema was observed with increasing AgNP concentration.

Our results show that compared to 2D cells in culture, which internalize extensively AgNW, the 3D RHE model showed no sign of penetration, which is consistent with the barrier function of the epidermis and the literature for many nanoparticles [START_REF] Crosera | Nanoparticle dermal absorption and toxicity: A review of the literature[END_REF][START_REF] Poland | Dermal Absorption of Nanomaterials; Danish Ministry of the Environment[END_REF]. Although, few examples of NP penetration into the viable layer of intact epidermis or even the dermis have been reported [START_REF] Leite-Silva | The effect of formulation on the penetration of coated and uncoated zinc oxide nanoparticles into the viable epidermis of human skin in vivo[END_REF][START_REF] Baroli | Penetration of metallic nanoparticles in human full-thickness skin[END_REF]. There is also strong evidence that the extent of nanomaterial penetration into the skin changes, when the integrity of the skin barrier is altered by several factors such as UV exposure, flexing, hydration status, skin age, chemical or physical damage, and skin pathologies [START_REF] Oberdörster | Nanotoxicology: An emerging discipline evolving from studies of ultrafine particles[END_REF][START_REF] Larese Filon | Occupational dermal exposure to nanoparticles and nano-enabled products: Part I-Factors affecting skin absorption[END_REF][START_REF] Schneider | Nanoparticles and their interactions with the dermal barrier[END_REF]. For example Mortensen et al. showed quantum dot nanoparticles penetration into the viable epidermis and dermis after exposure to UV simulating a sunburn [START_REF] Mortensen | In Vivo Skin Penetration of Quantum Dot Nanoparticles in the Murine Model: The Effect of UVR[END_REF]. Two studies found that even particles up to 1 µm in size can be taken up through intact skin to reach living cells under mechanical stress such as flexing [START_REF] Rouse | Effects of mechanical flexion on the penetration of fullerene amino acid-derivatized peptide nanoparticles through skin[END_REF][START_REF] Tinkle | Skin as a route of exposure and sensitization in chronic beryllium disease[END_REF]. Moreover, there are very few studies of dermal exposure to fibres or nanowires even though asbestos fibres have been shown to penetrate and persist in skin [START_REF] Langlands | Insulation workers in Belfast. 2. Morbidity in men still at work[END_REF].

Conclusions

The combination of 2D and 3D models for assessing the risk from dermal exposure to nanomaterials is expected to provide an increasingly important methodology as NW are incorporated into consumer products. Our results clearly show that keratinocytes, non-immune cell in the epidermis, rapidly and massively internalized AgNW leading to dose-dependent toxicity. In contrast, our study also demonstrates that the epidermis can provide an effective barrier for NW passage with no identified translocation or irritation end points even after a long exposure time (42 h). The present study, however, represents the first step in full risk assessment for acute dermal exposure to NW on intact epidermis and further studies of NW exposure using models for compromised skin and longer timescales with repeated exposure at low doses that are more realistic will provide a timely assessment of dermal exposure routes.

Supplementary Materials:

The following are available online at http://www.mdpi.com/2079-4991/8/4/232/s1. Electronic supplementary material (additional figures) as described in the text. Figure S1: SEM images showing the incomplete internalization of AgNW in human primary keratinocytes, Figure S2: SEM images of control non-treated primary keratinocytes, Figure S3: SEM images of cross-section of control non-treated primary keratinocytes, Table S1: Summary of the characteristics of two types of AgNW used in this study: short and long AgNW abbreviated S-AgNW and L-AgNW, Movie S1: 3D reconstruction of confocal images of primary human keratinocytes exposed to long AgNW for 24 h. Images were acquired in different mode: reflectance; transmittance; fluorescence in order to visualize respectively the NW; the cell and the nucleus that was stained with Hoechst. Scale bar represent 5 µm, Movie S2: Time-lapse confocal imaging of the kinetics of uptake of long AgNWs by primary human keratinocytes. Images were acquired each 2 min in reflected and transmitted light mode in order to visualize respectively the NW and the cells.

Figure 1 .

 1 Figure 1. Physicochemical behavior of AgNW in keratinocyte culture medium (KSFM). Kinetics of short and long AgNWs settling (A) and dissolution (B) following incubation in KSFM medium measured by UV-vis (A) and ICP-MS (B). (A) Series of optical absorption spectra acquired from unstirred 50 µg/mL AgNW suspensions in KSFM. Drop in absorbance at 380 nm reflects AgNW settling. (B) Ag + dissolution of a 50 µg/mL AgNW suspension over the time.

Figure 2 .

 2 Figure 2. (A) Viability of primary human keratinocytes cells following 24 h exposure with short, long AgNWs and AgNO 3 (represented each individually as bars from left to right respectively) determined by MTT assay. (B) Ion release control experiment. After 24 h incubation, cell medium containing AgNW or AgNO 3 were collected and centrifuged. Then the supernatants were collected and transferred onto other cells for 24 h prior to the MTT assay. Cell viability data is presented as normalized to the non-treated control cells (n = 3) +SD; Significant differences between non-treated and treated cells are indicated by a * where p < 0.05.

Figure 3 .

 3 Figure 3. Viability of primary human keratinocytes cells determined by MTT assay following 24 h treatment with short and long AgNWs (black and grey line respectively). Cell viability data is presented as normalized to the non-treated control cells (n = 3) + SD and represented using different dose metrics: mass of AgNW (A), number of NW (B), and NW surface area (C) relative to the cell growth surface area expressed in cm 2 .

Figure 4 .

 4 Figure 4. SEM images showing the uptake, the cellular distribution and changes of short (Left panel: A,B,E,F) and long (right panel: C,D,G,H) AgNWs in human primary keratinocytes cells after 24 h of exposure to 1.5 µg/mL AgNW. Images acquired at low (A-D) and high (E-H) magnification. SEM images were acquired using the secondary electron (SE; Left panel) and backscattered electron (BSE; Right panel) detector respectively. Scale bars represent 5 µm (A-D), 1 µm (E,F) and 0.4 µm (G,H) respectively. Changes in the morphology of AgNWs internalized inside cells is demonstrated (E,H) with bending (F) and a rough surface of internalized AgNW compared to their pristine surface outside cells (G,H). Insert in (H) show magnification of the box area and the modified AgNW surface.

Figure 5 .

 5 Figure 5. SEM images of cross sections of keratinocytes exposed to long AgNW for 72 h (A-D). Image B shows a higher magnification of the box area in (A). Sections of AgNW (highlighted by a white arrow) are localized in the cell cytoplasm, either free inside the cytoplasm or inside some endocytic vesicles. N = Nucleus, C = Cytoplasm, E/L = endosome/lysosome vesicles, ES = extracellular space. Scale bar represent 10 µm (A) and 2 µm (B-D).

Figure 6 .

 6 Figure 6. Confocal images showing the uptake and the cellular distribution of long AgNWs in human primary keratinocytes cells after 24 h of exposure to 1.5 µg/mL of L-AgNW. Scale bars represent 5 µm. Images were acquired in different mode (from left to right images): reflectance; transmittance; fluorescence in order to visualize respectively the NW; the cell and the nucleus that was stained with Hoechst. Image on the right is the merge of the three images. Images are shown at three different focal plans in order to visualize different depth in the cells.

Figure 7 .

 7 Figure 7. Hyperspectral analysis of AgNW and nanorods (AgNR) after 24 h of exposure in human primary keratinocytes (HPK). (A) Dark-field image of single HPK cell. (B) Region of HPK cell analyzed by visible-near infrared (V-NIR) hyperspectral microscopy. Circles indicate locations from where V-NIR spectra were extracted: white = reference spectrum; red = internal AgNW; green = internal AgNR. External AgNWs are outside this field of view (C) Processed V-NIR spectra compared with UV-vis absorbance from AgNW in HPK media.

Table 1 .Figure 8 .

 18 Figure 8. Phase contrast microscope (A-C) and TEM images (D,E) of reconstructed human epidermis (RHE) sections. No difference in RHE morphology was observed between non-treated (A) and s-AgNW (B) and L-AgNW (C) exposed RHE for 42 h. No AgNW was observed in the sections of RHE analysed (D,E), although the presence of AgNW lying at the surface of the stratum corneum or in the upper layer of the SC (E) is suspected (indicated by black arrow). The different layers of the skin are indicated. SC = stratum corneum; SG = stratum granulosum; SS = stratum spinosum; SB = stratum basale; M = membrane. The scale bar represent 20 µm (A-C); 10 µm (D) and 2 µm (E).

Acknowledgments:

The authors wish to thank the members of the Prométhée proteomic platform, M.J. Stasia and members of her team (IBP, CHU, Grenoble, France), K. Pernet Gallay (GIN, Grenoble, France), F. Roussel-Dherbey (CMTC, Grenoble, France), B. Faure and F. Saint Antonin (CEA, Grenoble, France). The NanoSetasizer and Ultramicrotome used in this study was part of Nano-ID platform which was funded by the EQUIPEX project ANR-10-EQPX-39-01. The project leading to this publication received funding from Excellence Initiative of Aix-Marseille University-A*MIDEX, a French "Investissements d'Avenir" program, through its associated Labex SERENADE project". This work was also supported from the Université Grenoble Alpes (UGA-PEPS AGIR) and the Labex OSUG@2020 (ANR10 LABX56).

Conflicts of Interest:

The authors declare no conflict of interest.

Abbreviations