1. Réduction du temps de recherche des lymphocytes NK en présence de spectateurs

Ce travail a été effectué au cours de mon stage de 5 mois du master 1 sous la direction de H. Rieger et K. Schwarz à Sarrebruck (Allemagne), pendant le premier semestre 2014. Il a donné lieu à une publication [P1].

Pendant ce stage, j'ai modélisé le processus de recherche des cellules tueuses naturelles (CTN) ou lymphocytes NK en présence de cellules spectatrices, en collaboration avec l'équipe de physique expérimentale de B. Qu à Hombourg (Allemagne). Les CTN cherchent à éliminer les cellules infectées par un pathogène (virus) ou tuméfiées, appelées ci-après cellules cibles, pour éviter la propagation de l'infection. Pendant leur recherche, elles rencontrent d'autres types de cellules non-cibles (saines) qui sont spectatrices mais qui augmentent l'efficacité et la migration de ces CTN. Cela est dû à la production de peroxyde d'hydrogène (H₂O₂) par ces cellules spectatrices. Les expériences in-vitro ont montré que la vitesse et la persistance des CTN étaient accrues en leur présence, sans destruction de leur capacité létale.

Nous avons alors développé des modèles diffusifs bidimensionnels, dans des domaines périodiques, permettant de comprendre les observations expérimentales. Les CTN ont été représentées comme des cellules browniennes discoïdes cherchant les cellules cibles immobiles instantanément détruites lorsqu'elles sont trouvées. Les cellules spectatrices, immobiles également, augmentent la diffusivité des CTN localement dans un rayon fini. Nous avons utilisé la méthode de Monte Carlo cinétique [1.1,1.2] d'une part, ainsi qu'un modèle sur réseau d'autre part pour reproduire cette dynamique stochastique numériquement. Ces simulations numériques montrent que le temps moyen pour localiser les cellules cibles, défini comme le temps de demi-vie des cibles, est réduit en présence des cellules spectatrices, autant par leur nombre que leur rayon d'action. Cet impact est supérieur à la présence de simples obstacles, même si ceux-ci favorisent déjà la recherche en diminuant la surface accessible. Notre étude numérique montre donc que la recherche des CTN est plus efficace en présence des cellules spectatrices, reproduisant parfaitement les observations des expériences in-vitro.

---

[P1] X. Zhou, R. Zhao, K. Schwarz, M. Mangeat, E. C. Schwarz, M. Hamed, I. Bogeski, V. Helms, H. Rieger, and B. Qu, Bystander cells enhance NK cytotoxic efficiency by reducing search time, Scientific Reports 7, 44357 (2017).[1.1] T. Oppelstrup, V. V. Bulatov, G. H. Gilmer, M. H. Kalos, and B. Sadigh, First-passage kinetic Monte Carlo method, Phys. Rev. Lett. 97, 230602 (2006).[1.2] K. Schwarz and H. Rieger, Efficient kinetic Monte Carlo method for reaction–diffusion problems with spatially varying annihilation rates, J. Comput. Phys. 237, 396 (2013).

2. Schémas d'approximation quantitatifs pour la transition vitreuse

Ce travail a été effectué au cours de mon stage de 2 mois du master 2 sous la direction de F. Zamponi au LPTENS de Paris, pendant l'hiver 2015. Il a donné lieu à ma première publication [P2].

Pendant ce stage, j'ai étudié les propriétés thermodynamiques de la transition vitreuse en développant un schéma d'approximation en partant de leurs solutions exactes en dimension infinie. La transition vitreuse est une transition du premier ordre aléatoire, appartenant à la classe d'universalité des verres de spin [2.1,2.2]. La phase vitreuse n'est donc pas unique, au contraire du cristal par exemple. Le diagramme de phase est alors extrêmement complexe, exhibant plusieurs états vitreux caractérisés par différentes propriétés thermodynamiques (par exemple l'équation d'état) qui dépendent généralement de la cinétique de la transition vitreuse.

Comme plusieurs modèles statistiques, les propriétés thermodynamiques de la transition vitreuse sont solvables en dimension infinie avec une structure de champ moyen [2.3]. Nous avons alors réalisé un développement systématique autour de cette solution en dimension infinie, similaire à un développement haute température ou basse densité. En utilisant la méthode des répliques [2.4], pour m copies du système original, nous avons alors dérivé des expressions analytiques pour l'entropie et l'équation d'état de la phase vitreuse en toutes dimensions d. Ces équations ne dépendent que de la thermodynamique et de la structure du liquide, et en particulier de la corrélation de paire g_liq et de l'entropie S_liq, et ont une structure comparable à celles obtenues avec la théorie du couplage de modes (TCM) [2.5] valable uniquement pour d=3. Seules des solutions approximatives de g_liq et S_liq existent dont les plus connues sont nommées Hyper-Netted-Chain et Percus-Yevick et dérivées à partir de l'équation du viriel [2.6].

À partir de l'équation d'état et de l'entropie (aussi appelée complexité Σ_eq) du verre, ainsi que des valeurs numériques de la corrélation de paire du liquide, nous avons calculé numériquement, en dimension d≥3, le facteur de non-ergodicité et les densités: (i) de la transition dynamique φ_d pour laquelle le liquide devient infiniment visqueux et plusieurs état vitreux apparaissent, et à partir de laquelle l'équation d'état pour m=1 admet des solutions; (ii) de la transition de Kauzmann φ_K pour laquelle le nombre d'état vitreux devient subexponentiel avec une crise entropique, Σ_eq=0 pour m=1, correspondant à la transition vitreuse idéale du second ordre; (iii) de la transition de jamming (m=0) allant de la densité de seuil φ_th à la densité de conditionnement compact du verre φ_GCP, correspondant respectivement aux transitions dynamique et de Kauzmann. Pour la dimension d=3, nous avons obtenu des valeurs numériques (φ_d≈0.53, φ_K≈0.62, φ_th≈0.45 et φ_GCP≈0.68) quasi-similaires aux précédentes études [2.7,2.8], basées notamment sur la TCM. Nous avons également obtenus des valeurs pour toutes dimensions d≥3. Ce travail a donc permis de relier les propriétés thermodynamiques de la transition vitreuse en dimension d=3 et en dimension infinie.

---

[P2] M. Mangeat and F. Zamponi, Quantitative approximation schemes for glasses, Phys. Rev. E 93, 012609 (2016).[2.1] D. Sherrington and S. Kirkpatrick, Solvable model of a spin-glass, Phys. Rev. Lett. 35, 1792 (1975).[2.2] T. R. Kirkpatrick, D. Thirumalai, and P. G. Wolynes, Scaling concepts for the dynamics of viscous liquids near an ideal glassy state, Phys. Rev. A 40, 1045 (1989).[2.3] G. Parisi and F. Zamponi, Mean-field theory of hard sphere glasses and jamming, Rev. Mod. Phys. 82, 789 (2010).[2.4] M. Mézard, G. Parisi, and M. Virasoro, Spin Glass Theory and Beyond, World Scientific (1987).[2.5] W. Götze, Complex Dynamics of Glass-Forming Liquids: A Mode-Coupling Theory, Oxford University Press, New York (2009).[2.6] J. P. Hansen and I. R. MacDonald, Theory of Simple Liquids, Academic Press, London (1986).[2.7] P. Charbonneau, Y. Jin, G. Parisi, and F. Zamponi, Hopping and the Stokes–Einstein relation breakdown in simple glass formers, PNAS 111, 15025 (2014).[2.8] G. Brambilla, D. El Masri, M. Pierno, L. Berthier, L. Cipelletti, G. Petekidis, and A. B. Schofield, Probing the equilibrium dynamics of colloidal hard spheres above the mode-coupling glass transition, Phys. Rev. Lett. 102, 085703 (2009).

3. Dispersion dans les milieux complexes

Ce travail a été effectué au cours de ma thèse réalisée sous la direction de D. S. Dean et T. Guérin au LOMA de Bordeaux, entre Octobre 2015 et Septembre 2018. Il a donné lieu à cinq publications [P3,P4,P5,P6,P7].

Pendant la première partie de ma thèse, j'ai analysé la dispersion de particules browniennes au sein de certains milieux complexes: les microcanaux et les milieux hétérogènes (réseaux d'obstacles). La diffusion dans les milieux simples, sans confinement entropique, est parfaitement caractérisée par les travaux réalisés au cours du XIX^e siècle posant les bases du mouvement brownien, des phénomènes diffusifs et de la mécanique statistique; et du début du XX^e siècle montrant que le mouvement brownien est un processus diffusif et posant les bases de la dynamique stochastique.

Or, la diffusion dans les milieux complexes est moins comprise, bien que largement étudiée. L'étude de la dispersion permet de comprendre l'étalement d'un nuage de particules browniennes initialement proches au cours du temps dans les milieux complexes. C'est un phénomène qui apparaît dans de nombreux contextes, comme le mélange ou le tri de particules, l'étalement de tâches de polluants, ou la cinétique des réactions chimiques. Nous nous sommes intéressés principalement à la dispersion dans les milieux confinés (microcanaux et réseaux d'obstacles) importante dans les systèmes biologiques [3.1], la microfluidique [3.2,3.3] ou les milieux poreux [3.4]. Pour ces systèmes, la distribution de Boltzmann n’est pas atteinte dans le régime des temps longs. Cependant, la taille typique de l'étalement des particules browniennes, donnée par le déplacement quadratique moyen, augmente au cours du temps et est caractérisée par la diffusivité effective $D_e$ inférieure à la diffusivité microscopique $D_0$, principal effet du confinement. Au cours de cette étude nous avons développé plusieurs expressions analytiques, valables dans différents régimes asymptotiques et vérifiées par des simulations numériques.

3.1 Dispersion et géométrie du confinement dans les microcanaux périodiques

Dans un premier temps, nous avons étudié la dispersion de particules browniennes dans les microcanaux périodiques, en nous intéressant uniquement au lien entre la géométrie du confinement et la dispersion. En effet, même pour des particules browniennes ponctuelles et sans interactions, dans un fluide au repos, il est bien connu que le confinement peut induire une réduction significative du ratio $D_e/D_0$. Les canaux considérés sont symétriques et bidimensionnels avec des parois situées en $y=\pm h(x)$ de période $L$, dans le système de coordonnées cartésiennes $(x,y)$. La diffusivité effective est alors reliée au déplacement quadratique moyen dans l'axe du canal $\langle x^2 \rangle \simeq 2 D_e t$ dans la limite des temps longs. La première théorie quantitative a été introduite par Jacobs [3.6], basée sur une réduction de dimensionnalité de l’espace et considérée aujourd'hui comme standard sous le nom de Fick-Jacobs (FJ). Elle permet de comprendre l'origine de cette diminution de la diffusivité. Une entropie effective $s(x) = k_B \ln h(x)$ peut être définie en comptant le nombre de configurations à la position $x$, en supposant un équilibre rapide dans la direction transverse. Les régions étroites et larges agissent alors respectivement comme des barrières et des pièges entropiques, ralentissant la dispersion des particules. La diffusivité effective est alors calculée en considérant les solutions exactes des processus diffusifs unidimensionnels [3.7]. De nombreuses études ont proposé des améliorations à cette approximation au cours des dernières décennies [3.8,3.9,3.10,3.11], en se basant pour la plupart sur une réduction de dimensionnalité introduisant une diffusivité unidimensionnelle $D(x)$, bien qu'elle mène généralement à des processus non-markoviens.

Tout d'abord, nous avons utilisé les formules de Kubo de T. Guérin et D. S. Dean [3.12,3.13], dérivées à partir de l'équation de Fokker-Planck et donnant une expression exacte de la diffusivité effective à partir d'une fonction auxiliaire $f(x,y)$, solution de l'équation de Laplace dans le canal. Ces formules permettent d'obtenir des résultats numériques exacts en déterminant la fonction auxiliaire par la méthode des éléments finis avec le logiciel FreeFem++ [3.5]; ainsi que des résultats analytiques asymptotiques. En considérant la fonction analytique $w(x+iy)= {\rm Re} f(x,y)$, nous avons montré que la diffusivité effective satisfait alors $${\rm Im} w\left(x+ih(x)\right) = h(x)-C, \qquad \frac{D_e}{D_0}= \frac{C}{\langle h \rangle}, \tag{3.1}$$ où la constante $C$ peut être déterminée grâce à la condition de périodicité. La notation $\langle h \rangle$ représente la moyenne spatiale de la hauteur $h(x)$ sur une période. Il s'agit d'un formalisme très compact, simplifiant les calculs.

Nous avons alors réalisé une classification des différents régimes de dispersion en fonction de la géométrie du canal. En définissant la hauteur minimale $a$ et la variation de hauteur $H$, la géométrie du canal est caractérisée par deux nombres sans dimensions : $\varepsilon = a/L$ et $\xi = H/a$ avec un profil donné par la hauteur $h(x) = a + H g(x/L)$ où $g(u)$ est une fonction $1$-périodique à valeurs dans l'intervalle $[0,1]$. Un exemple est donné dans l'encart de la Fig. 1a, pour lequel une trajectoire brownienne a été tracée. Dans la limite des canaux étroits ($\varepsilon \ll 1$), en définissant $h(x) = \varepsilon \zeta(x)$, nous avons dérivé à l'aide d'une méthode perturbative l'expression: $$\frac{D_e}{D_0} = \frac{1}{\langle h\rangle\langle h^{-1} \rangle } \left[1 - \frac{\varepsilon^2}{3} \frac{\langle \zeta'^2/\zeta\rangle}{\langle \zeta^{-1} \rangle} + {\cal O}(\varepsilon^4)\right]. \tag{3.2}$$ L'ordre dominant est donné par la diffusivité de FJ qui nous permet de distinguer deux types de canaux en fonction du comportement près de la hauteur minimale du canal, dans la limite de fortes ondulations $\xi \gg 1$. En supposant que la hauteur s'écrit $h(x) = a + \Lambda |x|^\nu$ près de son minimum, la diffusivité de FJ se réécrit en fonction du temps de premier passage à travers un entonnoir [3.14] pour un exposant $\nu>1$, et ne dépend pas de la géométrie du reste du canal. Pour un exposant $\nu < 1$, la diffusivité de FJ tends vers une constante dans la limite $\xi \rightarrow \infty$. Ces deux régimes de dispersion, contrôlés par la géométrie près de la hauteur minimale, définissent respectivement les canaux lisses et rugueux. De plus, nous avons montré que les résultats dérivés en considérant une réduction de dimensionnalité [3.8,3.9,3.10] sont bien valides, au moins jusqu'à l'ordre $\varepsilon^6$.

Dans la limite opposée des canaux larges ($\varepsilon \gg1$), qui n'avait jamais été abordée dans la littérature, la dispersion est dominée par les particules se situant dans la région $|y| < a$ où les particules contribuent au déplacement longitudinal alors que celles situées dans la région $|y| > a$ sont majoritairement piégées par les parois du canal. L'ordre dominant de la diffusivité effective est donc obtenu en considérant que les excursions dans les régions lentes n'y contribuent pas. Il s'agit du résultat observé pour la dispersion dans des peignes [3.15]. À partir de notre formalisme complexe et plusieurs transformations conformes du domaine périodique, nous avons dérivé l'ordre sous-dominant de la diffusivité effective qui ne dépend pas du profil du canal et faisant intervenir une constante universelle $\ln 2/\pi$: $$\frac{D_e}{D_0} = \frac{a}{\langle h\rangle} \left[ 1 + \frac{\ln 2}{\pi \varepsilon} + {\cal O}\left(\varepsilon^{-2}\right) \right]. \tag{3.3}$$

Dans les canaux de taille intermédiaire (pour des valeurs de $\varepsilon$ intermédiaires), nous avons défini un approximant de Padé permettant de relier les résultats asymptotiques pour $\varepsilon \ll1$ et $\varepsilon \gg 1$. Celui-ci approxime assez bien les valeurs numériques dans le régime intermédiaire pour les canaux lisses, mais présente quelques défauts pour les fortes ondulations ($\xi \gg1$). La dispersion est alors contrôlée par la présence de petites ouvertures. En considérant perturbativement les petites ouvertures [3.16,3.17], nous avons montré que la diffusivité effective s'écrit $$\frac{D_e}{D_0} \simeq \frac{\pi L}{2 \langle h \rangle} \frac{1}{\ln(2\kappa/\varepsilon)}, \tag{3.4}$$ où $\kappa$ dépend de la pseudo-fonction de Green [3.18] du domaine sans ouverture. Dans la limite $H \gg L$, cette constante a pour valeur $\kappa= 2/\pi$. À partir de ce résultat, la diffusivité effective peut s'écrire en fonction du temps de premier passage pour aller d'une petite ouverture, considérée réfléchissante, à la petite ouverture suivante [3.17]. Ce temps diffère ici du temps moyen de premier passage global pour atteindre une petite ouverture, par la correction du terme logarithmique donnée par $\kappa$. Nous avons alors obtenu trois régimes de dispersion contrôlés par la géométrie des canaux : canaux étroits, canaux larges et canaux présentant des petites ouvertures. Sur la Fig. 1a, nous avons représenté la diffusivité effective en fonction des paramètres $\varepsilon$ et $\xi$ pour le canal tracé en encart, en comparant les solutions numériques des formules de Kubo et les expressions analytiques obtenues pour ces trois régimes. Le diagramme de validité de ces expressions est montré sur la Fig. 1b, avec des lignes de transitions en $\varepsilon =1$ et $\varepsilon \xi^{1/\nu}=1$ pour les canaux lisses ($\nu ≥ 1$) ou $\varepsilon \xi=1$ pour les canaux rugueux ($\nu < 1$).

Enfin, nous avons abordé la dispersion dans les canaux étroits présentant des discontinuités de profil. Le résultat asymptotique développé pour les canaux étroits [Eq. (3.2)] n'est valable qu'à l'ordre dominant, dans la mesure où la dérivée de $\zeta(x)$ diverge. En utilisant notre formalisme complexe et en réalisant des transformations conformes du domaine, la diffusivité effective s'écrit pour $N$ discontinuités situées en $x=x_i$: $$\frac{D_e}{D_0} = \frac{1}{\langle h\rangle\langle h^{-1}\rangle} \left[ 1- \frac{\varepsilon}{L \langle \zeta^{-1} \rangle}\sum_{i=1}^N \gamma(\lambda_i) + {\cal O}(\varepsilon^2) \right], \tag{3.5}$$ où $\lambda_i$ est un paramètre relié à la $i^e$ discontinuité. Le piégeage des particules est donc plus important proche des discontinuités (terme d'ordre $\varepsilon$) que dans les sections continues (terme d'ordre $\varepsilon^2$), faisant diminuer tous deux la valeur de la diffusivité effective. La fonction $\gamma(\lambda)$ dépend de la géométrie de la discontinuité, et nous avons dérivé deux expressions exactes: $$\gamma_d(\lambda)=\frac{1+\lambda^2}{\pi\lambda} \ln \left\vert \frac{1+\lambda}{1-\lambda} \right\vert - \frac{2}{\pi} \ln \left\vert \frac{4\lambda}{1-\lambda^2}\right\vert; \qquad \gamma_c(\lambda)=-\frac{4}{\pi}\ln \left(\sin \frac{\pi \lambda}{2}\right), \tag{3.6}$$ où $\gamma_d(h_+/h_-)$ est valable pour une discontinuité simple avec une hauteur passant de $h_-$ à $h_+$ (canal 1 de la Fig. 1c), et $\gamma_c(h_0/h_-)$ pour des cloisons réduisant la hauteur au seul point de discontinuité $h_0 < h_-=h_+$ (canal 2 de la Fig. 1c). Sur la Fig. 1c, nous avons représenté la valeur de $\gamma$ obtenue à partir de la diffusivité effective, en comparant la solution numérique des formules de Kubo et des deux expressions analytiques dérivées. Nos résultats sont généralisables à d'autres formes de discontinuité. Il est à noter que ces fonctions $\gamma(\lambda)$ sont exactes et peuvent être reliées aux taux de piégeages introduits au niveau des discontinuités en considérant une réduction de dimensionnalité, calculés approximativement dans [3.19,3.20].

Nous avons donc décrit les régimes de dispersion pour toutes les géométries possibles d'un canal, et proposé une expression asymptotique de la diffusivité effective pour chacun d'eux en les reliant à celles utilisées dans la littérature. Ces expressions analytiques peuvent cependant être améliorées en calculant les ordres asymptotiques suivants.

3.2 Dispersion dans un réseau périodique d'obstacles attractifs

Dans un second temps, nous avons analysé la dispersion de particules browniennes ponctuelles dans un réseau périodique d'obstacles sphériques attractifs. La dispersion dans les milieux hétérogènes apparaît dans de nombreux contextes physiques, tels que les milieux poreux, les suspensions colloïdales ou les cellules vivantes, et permet de décrire les phénomènes de diffusion anormale [3.22]. La seule présence des obstacles sphériques (de rayon $R$) implique un ralentissement de la dispersion due à un piégeage entropique des particules. Il s’agit d’un mécanisme de piégeage entropique identique à celui mentionné dans la partie A. En dimension $d=2$, la dispersion dans un réseau d'obstacles est alors analogue à celle dans un canal, de même périodicité $L$ et de hauteur $h(x)=L/2 - \sqrt{R^2-x^2}$ [3.21], comme illustré sur la Fig. 2a. Cela nous permet de déterminer deux régimes spécifiques de dispersion dépendant uniquement de la fraction volumique des obstacles $\varphi = \pi R^2/L^2$. Dans la limite diluée ($\varphi \ll 1$), en utilisant le formalisme complexe développé pour les canaux [Eq. (3.1)], la diffusivité effective s'écrit $D_e/D_0 \simeq 1/(1+\varphi)$. Ce résultat est équivalent à celui de Maxwell, démontré pour la conductivité effective d’un milieu électrique [3.23,3.24]. Il est exact jusqu'à l'ordre $\varphi^4$ et est valide pour des fractions volumiques $\varphi < 0.4$. Dans la limite encombrée ($\varphi \simeq \pi/4$), l'approximation de FJ peut être appliquée. En remarquant que les formules de Kubo ne font intervenir que la normale à la sphère $n_x = x/R$, nous avons proposé une diffusivité unidimensionnelle $D(x) \simeq D_0 (1 - n_x^2/3)$ en remplaçant $h'(x)$ par $n_x$, équivalentes lors de la réduction de dimensionnalité [3.8,3.9,3.10]. Cette expression donne une meilleure approximation de la diffusivité effective que celles développées précédemment [3.21] et est valable pour des fractions volumiques à $\varphi > 0.67$. Sur la Fig. 2b, nous avons tracé la diffusivité effective en fonction de la fraction volumique $\varphi$, en comparant les solutions numériques des formules de Kubo (exactes) et les expressions analytiques mentionnées ci-dessus.

Ensuite, nous avons étudié l'effet d'un potentiel attractif $V(r)$ localisé sur la surface des sphères, sur la dispersion dans les réseaux d'obstacles. Les particules browniennes sont piégées dans les minima locaux de ce potentiel, entraînant un ralentissement bien connu et caractérisé par le temps moyen d'échappement du minimum de potentiel $t_{\rm esc} \sim \exp(\beta\Delta E)$, suivant la loi d'Arrhénius avec $\Delta E$ la barrière d'énergie à franchir et $\beta$ l'inverse de la température. Cependant, l'effet du couplage des deux barrières entropique (obstacles) et énergétique (potentiel attractif) sur la dispersion n'est pas si simple. Récemment, il a été mis en évidence par des simulations numériques [3.25] que la dispersion pouvait être accélérée pour un potentiel attractif, tout en respectant l'inégalité $D_e ≤D_0$. Ce résultat est contre-intuitif par le fait que la présence d'un minimum de potentiel puisse dépiéger une particule brownienne du confinement entropique des sphères. Cet effet pourrait cependant être compris si les surfaces des obstacles étaient suffisamment proches pour que les particules puissent passer de l'une à l'autre sans explorer tout l'espace. Mais ce phénomène apparaît également dans la limite diluée! Le potentiel que nous avons considéré est attractif $V(r)=-V_a <0$ sur $R<r<R_a$, présentant une barrière $V(r)=V_r>0$ sur $R_a<r<R_r$ et nul au delà (voir Fig. 2c). À l'aide des formules de Kubo [3.12,3.13], nous avons dérivé une expression de la diffusivité effective dans la limite diluée, cohérente avec celle obtenue pour un puits de potentiel carré ($V_r=0$) [3.26,3.27,3.28]. Nous avons alors remarqué que la dispersion était optimisée dans la limite de faible portée du potentiel $\alpha_a = (R_a-R)/R\ll 1$ et $\alpha_r = (R_r-R_a)/R\ll 1$, pour laquelle la diffusivité effective s'écrit en dimension $d=2$: $$\frac{D_e}{D_0} = 1- \frac{2(\tau+1)\xi^2 - (1-\tau)\xi +1}{(\tau+1)\xi+1} \varphi + {\cal O}(\varphi^2), \tag{3.7}$$ pour $\xi = \alpha_a \exp(\beta V_a)$ et $\tau = \alpha_a \alpha_r \exp[\beta (V_a+V_r)]$ fixés. Cette expression est tracée sur la Fig. 2d par rapport à la diffusivité effective maximale $D_e^*/D_0 \simeq 1 - 0.657 \varphi$. La dispersion est alors optimisée pour ce potentiel attractif si $D_e/D_0 > 1 - \varphi$, région définie par $\xi < 1/(1+\tau)$. De plus, les paramètres $\xi$ et $\tau$ peuvent être exprimés en fonction des temps de premier passage pour atteindre la surface $t_\Sigma$ en partant de l'extérieur (renormalisé par le volume) et pour s'échapper de celle-ci $t_{\rm esc}$ via les relations $\xi = t_{\rm esc}/t_{\Sigma}$ et $\tau = D_0 t_{\Sigma} / R^2$.

Enfin, nous avons étudié la dispersion des particules browniennes en considérant que la diffusion est médiée par la surface des sphères [3.29]. En considérant $D_\Sigma$ la diffusivité sur la surface, $D_{\rm B}$ la diffusivité dans le volume, $k$ le taux d'absorption sur la surface et $\lambda$ le taux de désorption, nous avons dérivé l'expression de la diffusivité effective à l'aide de formules de Kubo modifiées pour ce problème, par rapport aux équations de Fokker-Planck utilisées dans [3.29]. Cette expression revient à l'Eq. (3.7) dans la limite diluée pour $D_{\rm B}=D_\Sigma$, $\xi=k D_{\rm B}/\lambda R$ et $\tau=1/kR$. La dispersion dans un potentiel attractif localisé sur la surface des sphère est donc analogue à un processus de diffusion médiée par la surface des sphères et les taux d'absorption $k$ et de désorption $\lambda$ peuvent être reliés aux temps de premier passage définis précédemment: $t_{\rm esc} = \lambda^{-1}$ et ${t_\Sigma} = R/k D_{\rm B}$.

Nous avons donc décrit la dispersion dans un réseau d'obstacles en regardant l'effet entropique de ceux-ci d'une part, et le double effet entropique-énergique avec un potentiel attractif sur leur surface d'autre part. Le premier est parfaitement caractérisé tandis que le second a été analysé seulement dans la limite diluée montrant que la dispersion pouvait être accélérée par rapport à des obstacles réfléchissants.

---

[P3] M. Mangeat, T. Guérin, and D. S. Dean, Geometry controlled dispersion in periodic corrugated channels, EPL 118, 40004 (2017).[P4] M. Mangeat, T. Guérin, and D. S. Dean, Dispersion in two dimensional channels—the Fick–Jacobs approximation revisited, J. Stat. Mech. 2017, 123205 (2017).[P5] M. Mangeat, T. Guérin, and D. S. Dean, Dispersion in two-dimensional periodic channels with discontinuous profiles, J. Chem. Phys. 149, 124105 (2018).[P6] M. Mangeat, T. Guérin, and D. S. Dean, Effective diffusivity of Brownian particles in a two dimensional square lattice of hard disks, J. Chem. Phys. 152, 234109 (2020).[P7] A. Alexandre, M. Mangeat, T. Guérin, and D. S. Dean, How Stickiness Can Speed Up Diffusion in Confined Systems, Phys. Rev. Lett. 128, 210601 (2022).[3.1] P. C. Bressloff and J. M. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys. 85, 135–196 (2013).[3.2] X. Yang, C. Liu, Y. Li, F. Marchesoni, P. Hänggi, and H. P. Zhang, Hydrodynamic and entropic effects on colloidal diffusion in corrugated channels, PNAS 114, 9564–9569 (2017).[3.3] M. Aminian, F. Bernardi, R. Camassa, D. M. Harris, and R. M. McLaughlin, How boundaries shape chemical delivery in microfluidics, Science 354, 1252 (2016).[3.4] H. Brenner, Dispersion resulting from flow through spatially periodic porous media, Phil. Trans. R. Soc. Lond. A 297, 81–133 (1980).[3.5] F. Hecht, New development in freefem++, J. Num. Math. 20, 251–266 (2013).[3.6] M. H. Jacobs, Diffusion Processes, Springer, Berlin, Heidelberg (1935).[3.7] S. Lifson and J. L. Jackson, On the Self-Diffusion of Ions in a Polyelectrolyte Solution, J. Chem. Phys. 36, 2410 (1962).[3.8] R. Zwanzig, Diffusion Past an Entropy Barrier, J. Phys. Chem. 96, 3926 (1991).[3.9] D. Reguera and J. M. Rubí, Kinetic equations for diffusion in the presence of entropic barriers, Phys. Rev. E 64, 061106 (2001).[3.10] P. Kalinay and J. K. Percus, Corrections to the Fick-Jacobs equation, Phys. Rev. E 74, 041203 (2006).[3.11] K. D. Dorfman and E. Yariv, Assessing corrections to the Fick–Jacobs equation, J. Chem. Phys. 141, 044118 (2014).[3.12] T. Guérin and D. S. Dean, Force-Induced Dispersion in Heterogeneous Media, Phys. Rev. Lett. 115, 020601 (2015).[3.13] T. Guérin and D. S. Dean, Kubo formulas for dispersion in heterogeneous periodic nonequilibrium systems, Phys. Rev. E 92, 062103 (2015).[3.14] D. Holcman and Z. Schuss, Control of flux by narrow passages and hidden targets in cellular biology, Rep. Prog. Phys. 76, 074601 (2013).[3.15] D. S. Dean and K. M. Jansons, Brownian excursions on combs, J. Stat. Phys. 70, 1313–1332 (1993).[3.16] M. J. Ward and J. B. Keller, Strong localized perturbations of eigenvalue problem, SIAM J. Appl. Math. 53, 770 (1993).[3.17] S. Pillay, M. Ward, A. Peirce, and T. Kolokolnikov, An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part I: Two-Dimensional Domains, Multiscale Model. Simul. 8, 803 (2010).[3.18] G. Barton, Elements of Green’s Functions and Propagation: Potentials, Diffusion, and Waves, Clarendon Press (1989).[3.19] A. M. Berezhkovskii, M. I. Monine, C. B. Muratov, and S. Y. Shvartsman, Homogenization of boundary conditions for surfaces with regular arrays of traps, J. Chem. Phys. 124, 036103 (2006).[3.20] P. Kalinay and J. K. Percus, Mapping of diffusion in a channel with abrupt change of diameter, Phys. Rev. E 82, 031143 (2010).[3.21] L. Dagdug, M.-V. Vazquez, A. M. Berezhkovskii, V. Yu. Zitserman, and S. M. Bezrukov, Diffusion in the presence of cylindrical obstacles arranged in a square lattice analyzed with generalized Fick-Jacobs equation, J. Chem. Phys. 136, 204106 (2012).[3.22] M. J. Saxton, Anomalous Diffusion Due to Obstacles: A Monte Carlo Study, Biophys. J 92, 1178–1191 (2007).[3.23] J. C. Maxwell, A Treatise on Electricity and Magnetism, Clarendon Press (1881).[3.24] J. R. Kalnin, E. A. Kotomin, and J. Maier, Calculations of the effective diffusion coefficient for inhomogeneous media, J. Phys. Chem. Solids 63, 449 (2001).[3.25] G. G. Putzel, M. Tagliazucchi, and I. Szleifer, Nonmonotonic Diffusion of Particles Among Larger Attractive Crowding Spheres, Phys. Rev. Lett. 113, 138302 (2014).[3.26] B. Cichocki and B. U. Felderhof, Diffusion coefficients and effective viscosity of suspensions of sticky hard spheres with hydrodynamic interactions, J. Chem. Phys. 93, 4427 (1990).[3.27] B. Cichocki and B. U. Felderhof, Self-diffusion in suspensions of interacting Brownian particles, Phys. Rev. A 42, 6024 (1990).[3.28] B. Cichocki and B. U. Felderhof, Time-Dependent Self-Diffusion of Brownian Particles with Square Well Interaction, Langmuir 8, 2889 (1992).[3.29] J.-F. Rupprecht, O. Bénichou, D. S. Grebenkov, and R. Voituriez, Exact mean exit time for surface-mediated diffusion, Phys. Rev. E 86, 041135 (2012).

4. Vortex browniens créés par les pièges optiques

Ce travail a été effectué au cours de ma thèse réalisée sous la direction de D. S. Dean et T. Guérin au LOMA de Bordeaux, entre Octobre 2015 et Septembre 2018. Il a donné lieu à trois publications [P8,P9,P10].

Pendant la seconde partie de ma thèse, en collaboration avec Y. Louyer et Y. Amarouchene pour la partie expérimentale, j'ai analysé la dynamique stochastique de particules browniennes dans un piège optique en présence d'une force non conservative. Ces particules sont soumises à un piège purement énergétique, contrairement à celles de la partie 3 où elles étaient confinées par un piège entropique. Les pièges (ou pinces) optiques développés par Ashkin dans les années 1970 [4.2] permettent de manipuler des objets de petites dimensions (cellules vivantes, colloïdes, nanoparticules) à l'aide d'un faisceau laser et de mesurer l'intensité de petites forces. Ils ont permis de caractériser l'élasticité des brins d'ADN et vérifier des résultats de la physique statistique hors-équilibre et de la mécanique quantique. Le piège optique y est essentiellement supposé harmonique, or celui-ci est généralement anharmonique et la pression de radiation du laser crée une force non conservative dans sa direction longitudinale [4.3]. Le système est donc hors-équilibre et l'état stationnaire ne vérifie pas la distribution de Boltzmann, conduisant à la présence de courants. Les particules suivent alors une circulation toroïdale dans ce flux stationnaire, appelée vortex brownien [4.4] et observée à températures et pressions ambiantes, caractéristiques du régime suramorti [4.5].

Nous nous sommes intéressés au régime sous-amorti qui n'avait jamais été étudié, avec un coefficient de friction $\gamma$ faible pour les pressions atteintes dans les expériences menées en parallèle (i.e. $P<10 {\rm mbar}$). L'équation de Langevin sous-amortie pour la trajectoire ${\bf x}(t)$ de la particule brownienne en dimension $d=3$ est $$m \ddot {\bf x} = -\gamma \dot {\bf x} +{\bf F}({\bf x}) +\sqrt{2 k_B T\gamma}{\bf \xi}, \tag{4.1}$$ avec $m$ la masse de la particule brownienne et $T$ la température. ${\bf \xi}$ représente le bruit thermique, considéré comme un bruit blanc gaussien. La force ${\bf F}({\bf x})$ créée par le piège optique [4.3] est composée d'une partie conservative dérivant du potentiel harmonique $V_{\rm harm}({\bf x}) = \kappa \left(x^2+ y^2 + \eta z^2\right)/2$ avec une raideur $\kappa$ dans les directions transverses et $\eta \kappa$ dans la direction longitudinale ($\eta \simeq 0.2$); et d'une partie non-conservative ${\bf F}_{\rm scat}({\bf x}) = \varepsilon \kappa a [ 1-(x^2+y^2)/a^2 ] {\bf e_z}$, appelée force de scattering, pour une largeur de faisceau $a$ et d'amplitude $\varepsilon$. Nous avons réalisé des simulations numériques [4.1] en intégrant l'équation de Langevin sous-amortie [Eq. (4.1)], permettant d'obtenir l'état stationnaire des vortex browniens pour toute amplitude $\varepsilon$ et tout facteur de qualité $Q= \sqrt{\kappa m}/\gamma$.

Tout d'abord, nous avons dérivé la densité de probabilité stationnaire dans le régime sous-amorti en considérant la force de scattering comme une perturbation ($\varepsilon \ll 1$), et donc une déviation à la distribution de Maxwell-Boltzmann $P_0({\bf x},{\bf v})$ dans l'espace des phases $({\bf x},{\bf v})$. En intégrant sur les vitesses ${\bf v}$, nous avons alors obtenu la densité de probabilité marginale en positions ${\bf x}$ et l'expression du courant, formant le vortex brownien, dans l'espace des positions ${\bf x}$: $$\frac{{\bf J} ({\bf x})}{P_0({\bf x})} \simeq \varepsilon \omega_0 A_x(Q) \sqrt{8\pi \beta \kappa} \left[ \eta z (x {\bf e_x} + y {\bf e_y} )+\left(\frac{2}{\beta\kappa}- x^2-y^2\right) {\bf e_z} \right], \tag{4.2}$$ avec $\omega_0 = \sqrt{\kappa/m}$ la pulsation caractéristique du piège harmonique. La géométrie des vortex est indépendante du taux d'amortissement caractérisé par $Q$. Sur la Fig. 3a, nous avons représenté la densité de particule et la géométrie du vortex pour $\varepsilon = 0.1$ et $Q=1$. Leur amplitude $A_x(Q)$, proportionnelle à la circulation du courant $\Omega_0$, présente un maximum en fonction du facteur de qualité et est tracée sur la Fig. 3c pour $\varepsilon=0.1$. Le courant dans l'espace des vitesses se présente également sous la forme d'un vortex dont seule l'amplitude dépend de $Q$, avec une évolution strictement croissante.

Ensuite, nous avons dérivé l'expression de la densité spectrale $S_{zz}(\omega)$ correspondant à la transformée de Fourier de la corrélation $\langle z(t) z(0) \rangle$. Le processus $z$ peut être décomposé en deux processus indépendants : $z_{\rm B}$ soumis uniquement au bruit thermique et $z_{\rm nc}$ soumis à la force non conservative faisant intervenir les carrés des processus $x$ et $y$, dont les densités spectrales ont une forme lorentzienne avec un maximum en $\omega=\omega_0$ et une largeur $\Gamma = \gamma/m$ à demi-hauteur. Cette dernière est tracée sur la Fig. 3d pour un facteur de qualité $Q=100$. Le processus $z_{\rm B}$ à une densité spectrale similaire. Après avoir calculé la fonction d'autocorrélation de $x^2$, l'expression de la densité spectrale du processus $z$ s'écrit $$S_{zz}(\omega)= \frac{k_B T}{m} \frac{2 \Gamma}{(\eta\omega_0^2 -\omega^2)^2 + \Gamma^2 \omega^2} \left[1 + \frac{8\varepsilon^2 \omega_0^2}{a^2} \frac{k_B T}{m} \frac{4\Gamma^2 + 4 \omega_0^2 + \omega^2}{(\Gamma^2 + \omega^2)[(4\omega_0^2- \omega^2)^2+4 \Gamma^2 \omega^2]}\right]. \tag{4.3}$$ Deux pics sont présents en $\omega= \sqrt \eta \omega_0$ et $\omega = 2 \omega_0$, correspondant à la signature respective des processus conservatif et non-conservatif. La force non-conservative influe également sur l'état stationnaire via le plateau situé en $\omega=0$, se comportant comme $1/\gamma$ dans la limite sous-amortie. Cette densité spectrale est également tracée sur la Fig. 3d, pour le facteur de qualité $Q=100$.

Pour finir, nous avons dérivé une expression exacte de la densité de probabilité stationnaire $P_s({\bf x})$ dans le régime suramorti, qui n'avait été calculé que dans la limite $\varepsilon \ll 1$ [4.5]. Les processus d'Ornstein-Uhlenbeck $x$, $y$ et $z_{\rm B}$ sont décrits par la distribution de Boltzmann dans l'état stationnaire. Seul le processus $z_{\rm nc}$ est hors-équilibre. Sa densité de probabilité $P_{\rm nc}$ peut être écrite sous la forme d'une intégrale de chemin sur les trajectoires gaussiennes $x(t)$ et $y(t)$. En utilisant les propriétés du processus d'Ornstein-Uhlenbeck, $P_{\rm nc}$ s'écrit dans l'espace de Laplace de variable $s$ comme le déterminant d'un opérateur que nous avons pu calculer [4.6,4.7]. Nous avons alors montré que $\widehat P_{\rm nc}(x,y,s)$ à une forme gaussienne en $x$ et $y$, de norme $\widehat P_{\rm nc}(s)$ et de variance $v(s)$, tels que $$\widehat P_{\rm nc}(s) = \frac{(X_s/2)^{\alpha -1 }}{\Gamma(\alpha) I_{\alpha-1}(X_s)} \qquad {\rm et} \qquad v(s) = \frac{2\alpha}{(\beta \kappa)^2}\frac{I_\alpha(X_s)}{X_s I_{\alpha-1} (X_s)}, \tag{4.4}$$ avec $X_s = (4/\eta\beta\kappa) \sqrt{\varepsilon s/a}$ et $\alpha=2/\eta$. Dans l'espace réel, la densité de probabilité stationnaire $P_s({\bf x})$ s'écrit comme la convolution de $\widehat P_{\rm nc}$ avec $\widehat P_{\rm B}$, la transformée de Fourier de la distribution de Boltzmann du processus $z_{\rm B}$. Cette expression exacte nous permet de retrouver la théorie de perturbation d'ordre $\varepsilon$ et d'obtenir l'ordre suivant, autant pour la densité de probabilité que pour le courant. Sur la Fig. 3b, nous avons tracé la densité de particule et la géométrie du vortex pour $\varepsilon = 10$ et $Q=0.1$, obtenues avec les simulations numériques et identiques aux résultats analytiques.

Nous avons donc caractérisé les vortex browniens dans le régime sous-amorti pour des faibles amplitudes et dans le régime suramorti pour tout valeur d'amplitude, en dérivant une expression analytique pour la densité de probabilité, le courant et la densité spectrale en accord avec les observations expérimentales. Seul le régime sous-amorti pour des grandes amplitudes n'a pas été complètement analysé.

---

[P8] Y. Amarouchene, M. Mangeat, B. Vidal Montes, L. Ondic, T. Guérin, D. S. Dean, and Y. Louyer, Nonequilibrium Dynamics Induced by Scattering Forces for Optically Trapped Nanoparticles in Strongly Inertial Regimes, Phys. Rev. Lett. 122, 183901 (2019).[P9] M. Mangeat, Y. Amarouchene, Y. Louyer, T. Guérin, and D. S. Dean, Role of nonconservative scattering forces and damping on Brownian particles in optical traps, Phys. Rev. E 99, 052107 (2019).[P10] M. Mangeat, T. Guérin, and D. S. Dean, Steady state of overdamped particles in the non-conservative force field of a simple non-linear model of optical trap, J. Stat. Mech. 2021, 113205 (2021).[4.1] D. A. Sivak, J. D. Chodera, and G. E. Crooks, Time Step Rescaling Recovers Continuous-Time Dynamical Properties for Discrete-Time Langevin Integration of Nonequilibrium Systems, J. Phys. Chem. B 118, 6466–6474 (2014).[4.2] A. Ashkin, Acceleration and Trapping of Particles by Radiation Pressure, Phys. Rev. Lett. 24, 156–159 (1970).[4.3] J. Gieseler, L. Novotny, and R. Quidant, Thermal nonlinearities in a nanomechanical oscillator, Nat. Phys. 9, 806 (2013).[4.4] B. Sun, J. Lin, E. Darby, A. Y. Grosberg, and D. G. Grier, Brownian vortexes, Phys. Rev. E 80, 010401 (2009).[4.5] H. W. Moyses, R. O. Bauer, A. Y. Grosberg, and D. G. Grier, Perturbative theory for Brownian vortexes, Phys. Rev. E 91, 062144 (2015).[4.6] K. Kirsten and A. J. McKane, Functional determinants by contour integration methods, Annals of Physics, Annals of Physics 308, 502–527 (2003).[4.7] K. Kirsten and A. J. McKane, Functional determinants for general Sturm–Liouville problems, J. Phys. A: Math. Gen. 37, 4649 (2004).

5. Temps d'échappement d'une cellule et organisation du cytosquelette

Ce travail a été effectué au cours de mon post-doctorat dans le groupe de H. Rieger à Sarrebruck (Allemagne), commencé en Novembre 2018. Il a donné lieu à deux publications [P11,P12].

Au cours de mon post-doctorat, j'ai analysé le temps mis par une particule brownienne pour s'échapper d'une cellule vivante, en considérant un modèle simplifié de l'organisation de son cytosquelette. Le transport intracellulaire de diverses particules-cargos (protéines, vésicules, mitochondries par exemple) vers la synapse immunologique située sur la membrane cellulaire, est crucial pour le bon fonctionnent des cellules et des organismes. Le cytosquelette est auto-organisé par un réseau de filaments, majoritairement constitué par des microtubules polarisés et des filaments d'actine, transportant rapidement les particules entre les différentes régions de la cellule. Ces deux types de filaments séparent le cytosquelette en deux régions : les microtubules émergeant radialement du centrosome, et les filaments d'actine distribués aléatoirement formant un fin cortex près de la membrane cytoplasmique [5.1]. Les particules-cargos réalisent une recherche intermittente [5.2,5.3], en alternant aléatoirement entre un transport diffusif dans le cytoplasme et un transport balistique sur les filaments [5.4]. Ce transport balistique est assisté par des moteurs moléculaires liant les particules-cargos aux filaments: la dynéine (du centrosome au cortex) et la kinésine (du cortex au centrosome) sur les microtubules, et la myosine sur les filaments d'actine. Il a été montré récemment que le temps pour atteindre la synapse via cette recherche intermittente pouvait être minimisé en jouant sur la taille du cortex, grâce à une diffusion effective plus importante dans celui-ci [5.5,5.6].

Nous nous sommes alors intéressés au calcul du temps moyen de premier passage (MFPT pour mean first passage time) pour atteindre la synapse immunologique en considérant une géométrie simplifiée du cytosquelette, représentée sur la Fig. 4a, pour des particules passives. Ces particules réalisent un mouvement diffusif dans un disque de rayon $R=1$ composé de deux régions concentriques, avec une diffusivité $D_2$ dans la région centrale et $D_1>D_2$ dans le cortex de taille $\Delta$, cohérent avec les études récentes [5.5,5.6]. De plus, une différence de potentiel $\Delta V = V_2 - V_1$ est présente entre ces deux régions, et la synapse est modélisée par un arc de cercle absorbant de taille $\varepsilon$. En partant d'une coordonnée ${\bf x}$, le MFPT est noté $t_1({\bf x})$ dans le cortex et $t_2({\bf x})$ dans la région centrale. Il satisfait une équation de Poisson dans chacune des deux régions: $D_i \nabla^2 t_i = -1$, avec une condition aux limites mixte: $t_1= 0$ au niveau de la synapse absorbante et ${\bf n} \cdot \nabla t_1= 0$ sur le reste de la membrane réfléchissante, avec ${\bf n}$ la normale à la surface, et avec les conditions: $t_1 = t_2$ et $D_1 e^{-\beta V_1}{\bf n} \cdot \nabla t_1 = D_2 e^{-\beta V_2}{\bf n} \cdot \nabla t_2$ à l'interface des deux régions [5.7]. Depuis plusieurs décennies, le MPFT a été largement analysé pour des petites régions d'échappement [5.8,5.9,5.10] correspondant ici à la limite $\varepsilon \ll 1$. Dans le cas d'un disque homogène ($D_1=D_2$ et $\Delta V=0$), une expression du MFPT a été dérivée dans cette limite de petite ouverture [5.11] ainsi qu'une expression exacte [5.12] démontrée à partir de transformations conformes. Dans la suite, je ne parlerai que des résultats obtenus en 2d, tous reproductibles en 3d et présentant qualitativement les mêmes conclusions.

Nous avons dérivé une expression analytique du MFPT dans la limite $\varepsilon \ll 1$ en supposant la petite ouverture comme une perturbation [5.8,5.9]. En particulier, nous avons étudié le temps moyen de premier passage global (GMFPT): $T = D_1 \langle t \rangle/R^2$ pour lequel la position de départ est moyennée par rapport à la densité de probabilité stationnaire, satisfaisant la distribution de Boltzmann. Il s'écrit sous la forme $$\displaylines{T \simeq \left[1-\chi^2(1-\xi)\right] \left[- \ln \frac{\varepsilon}{4} - 2 \sum_{k=1}^\infty \left( \frac{D_1-D_2\xi}{D_1 + D_2 \xi} \right)^k \ln(1-\chi^{2k}) \right] + \frac{1}{8} \\+ \frac{(D_1-D_2)\chi^2 \xi - 3D_2 (1-\xi)(1-\chi^2)}{8D_2 \left[1-\chi^2(1-\xi)\right]} \chi^2 - \frac{(1-\xi)^2 \chi^4\ln \chi}{2\left[1-\chi^2(1-\xi)\right]},} \tag{5.1}$$ avec $\chi = 1 - \Delta$ et $\xi= \exp(-\beta \Delta V)$. De même, nous avons dérivé une expression du MFPT dans la limite d'un fin cortex ($\Delta \ll 1$), similaire au problème de diffusion médiée par la surface [5.13]. Dans cette limite, le GMFPT s'écrit $$ T \simeq \begin{cases} \frac{(2\pi-\varepsilon)^3}{24\pi}, \qquad &\Delta V = +\infty\\ \frac{D_1}{D_2} \left( - \ln \frac{\varepsilon}{4} + \frac{1}{8} \right), \qquad &\Delta V < +\infty \end{cases}\tag{5.2}$$ où $\Delta V = + \infty$ signifie que la région centrale est inaccessible, avec une interaction de cœur dur.

Lorsque les particules ne peuvent pas entrer dans la région centrale ($\Delta V= +\infty$), nous observons que ce temps global peut être minimisé pour un cortex de taille $\Delta \ne 1$. Sur la Fig. 4b, nous avons tracé le GMFPT calculé numériquement avec FreeFem++ [5.14] à partir des équations aux dérivées partielles du MFPT, donnant des résultats similaires aux simulations numériques réalisées à l'aide de la méthode de Monte Carlo cinétique [5.15,5.16]. L'expression donnée par l'Eq. (5.1) est valable pour des cortex pas trop petits, définis par $-\Delta \ln \varepsilon>1$, englobant le minimum. Le raccord des expressions du GMFPT dans les limites $\varepsilon \ll 1$ et $\Delta \ll1$ donne la valeur asymptotique $T \simeq - 2 \Delta \ln \varepsilon + \pi^2/3$.

Sans différence de potentiel ($\Delta V=0$), nous avons montré que le GMFPT était une fonction strictement décroissante de l'angle $\varepsilon$ et de la taille du cortex $\Delta$ et strictement croissante du ratio des diffusivités $D_1/D_2$. Aucune optimisation ne peut alors être réalisée. Cela peut s'expliquer par l'absence de biais forçant les particules à aller vers le cortex, comme c'était le cas avec un cœur répulsif et lors des études menées avec une diffusion médiée par la surface [5.13] ou avec un transport balistique intermittent [5.5,5.6]. Malgré cela, nous avons pu montré une propriété intéressante due à l'hétérogénéité de la diffusivité. Le point de départ pour lequel le MFPT est maximal dans le disque évolue discontinûment avec le ratio $D_1/D_2$. Sa distance $r_{\rm max}$ avec le centre du domaine est tracée sur la Fig. 4c. Lorsque $D_1 = D_2$, il est évident que le point opposé à l'ouverture ($r_{\rm max}=1$) possède le temps maximal pour atteindre l'ouverture. Or, la région centrale devient une région lente en augmentant la diffusivité du cortex et le centre devient le point avec un MFPT maximal ($r_{\rm max}=0$) dans la limite $D_1 \gg D_2$. La transition est alors discontinue avec $r_{\rm max}$ passant de $1$ à une valeur inférieure à $1-\Delta$ pour un ratio de diffusivité déterminé sur la Fig. 4c par la courbe en pointillé.

Dans le cas général, nous avons montré que GMPFT présentait un minimum pour $\Delta \ne 1$ si et seulement si la différence de potentiel satisfait l'inégalité (dans la limite $\varepsilon \ll 1$): $$\beta \Delta V \gtrsim \frac{2(D_1-D_2)}{(D_1+D_2)\left(-\ln \frac{\varepsilon}{4}+\frac{3}{8}\right)}.\tag{5.3}$$ Sur la Fig. 4d, nous avons tracé le GMFPT calculé numériquement avec FreeFem++ [5.14] à partir des équations aux dérivées partielles du MFPT pour différentes différences de potentiel. Nous observons qu'un minimum est présent pour $\Delta \ne 1$ seulement si $\beta \Delta V \gtrsim 0.5$. L'expression donnée par l'Eq. (5.1) est valable pour des cortex pas trop petits, englobant ces minima. Le GMFPT tend vers une constante indépendante de la différence de potentiel dans la limite $\Delta \ll 1$, donnée par l'Eq. (5.2). Sur la Fig. 4e, nous avons représenté le diagramme d'optimisation du GMFPT déduit à partir des solutions numériques obtenues avec FreeFem++ et de la condition analytique donnée par l'Eq. (5.3). Nous pouvons remarqué que pour $\beta \Delta V > 0.68$, le GMFPT présente un minimum quel que soit le ratio des diffusivités $D_1/D_2$.

Nous avons donc dérivé des expressions analytiques pour le MFPT pour une diffusivité hétérogène dans la limite des petites ouvertures, et montré qu'il était nécessaire d'avoir un mécanisme forçant les particules à aller vers la membrane cellulaire pour minimiser le MFPT avec la taille du cortex d'actine.

---

[P11] M. Mangeat and H. Rieger, The narrow escape problem in a circular domain with radial piecewise constant diffusivity, J. Phys. A: Math. Theor. 52, 424002 (2019).[P12] M. Mangeat and H. Rieger, Narrow escape problem in two-shell spherical domains, Phys. Rev. E 104, 044124 (2021).[5.1] B. Alberts, A. Johnson, J. Lewis, D. Morgan, M. Raff, K. Roberts, and P. Walter, Molecular Biology of the Cell, Garland Science, New York (2014).[5.2] C. Loverdo, O. Benichou, M. Moreau, and R. Voituriez, Enhanced reaction kinetics in biological cells, Nat. Phys. 4, 134-137 (2008).[5.3] O. Bénichou, C. Loverdo, M. Moreau, and R. Voituriez, Intermittent search strategies, Rev. Mod. Phys. 83, 81-129 (2011).[5.4] P. C. Bressloff and J. M. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys. 85, 135–196 (2013).[5.5] K. Schwarz, Y. Schröder, B. Qu, M. Hoth, and H. Rieger, Optimality of Spatially Inhomogeneous Search Strategies, Phys. Rev. Lett. 117, 068101 (2016).[5.6] A. E. Hafner and H. Rieger, Spatial organization of the cytoskeleton enhances cargo delivery to specific target areas on the plasma membrane of spherical cells, Phys. Biol. 13, 066003 (2016).[5.7] S. Redner, A Guide to First-Passage Processes, Cambridge University Press (2001).[5.8] M. J. Ward and J. B. Keller, Strong localized perturbations of eigenvalue problem, SIAM J. Appl. Math. 53, 770 (1993).[5.9] S. Pillay, M. Ward, A. Peirce, and T. Kolokolnikov, An Asymptotic Analysis of the Mean First Passage Time for Narrow Escape Problems: Part I: Two-Dimensional Domains, Multiscale Model. Simul. 8, 803 (2010).[5.10] O. Bénichou and R. Voituriez, Narrow-Escape Time Problem: Time Needed for a Particle to Exit a Confining Domain through a Small Window, Phys. Rev. Lett. 100, 168105 (2008).[5.11] A. Singer, Z. Schuss, and D. Holcman, Narrow Escape, Part II: The Circular Disk, J. Stat. Phys. 122, 465 (2006).[5.12] C. Caginalp and X. Chen, Analytical and Numerical Results for an Escape Problem, Arch Rational Mech Anal 203, 329–342 (2012).[5.13] J.-F. Rupprecht, O. Bénichou, D. S. Grebenkov, and R. Voituriez, Exact mean exit time for surface-mediated diffusion, Phys. Rev. E 86, 041135 (2012).[5.14] F. Hecht, New development in freefem++, J. Num. Math. 20, 251–266 (2013).[5.15] T. Oppelstrup, V. V. Bulatov, G. H. Gilmer, M. H. Kalos, and B. Sadigh, First-passage kinetic Monte Carlo method, Phys. Rev. Lett. 97, 230602 (2006).[5.16] K. Schwarz and H. Rieger, Efficient kinetic Monte Carlo method for reaction–diffusion problems with spatially varying annihilation rates, J. Comput. Phys. 237, 396 (2013).

6. Mouvement collectif sur réseaux

Ce travail a été effectué au cours de mon post-doctorat dans le groupe de H. Rieger à Sarrebruck (Allemagne), commencé en Novembre 2018. Il a donné lieu à quatre publications [P13,P14,P15,P16].

Pendant mon post-doctorat, en collaboration avec R. Paul et S. Chatterjee de l'IACS de Calcutta (Inde), j'ai également étudié le mouvement de particules actives s'alignant localement avec une interaction de Potts. Les mouvements collectifs sont largement observés dans la nature (nuées d'oiseaux, troupeaux de mammifères, bancs de poissons, essaims de bactéries) et également étudiés dans divers systèmes artificiels (gouttes de liquides, colloïdes roulants, disques polarisés vibrants) [6.1]. Ils sont constitués d'un grand nombre de particules actives qui consomment de l'énergie pour s'autopropulser et interagissent en s'alignant entre elles. Un mouvement synchrone de grands amas de particules, appelés nuées, émerge spontanément pour des grandes densités ou des faibles fluctuations thermiques. Cette transition vers les nuées (flocking transition en anglais) est un phénomène hors-équilibre décrit par de nombreux modèles théoriques. Le premier d'entre eux est le modèle de Vicsek (MV) [6.2,6.3] pour lequel les particules actives réalisent un mouvement balistique avec une vitesse constante dans une direction s'alignant avec celles de leurs voisines. Cet alignement est similaire à une interaction ferromagnétique dans un bain thermique, et un mouvement collectif émerge grâce à la brisure spontanée de la symétrie de rotation. Une théorie hydrodynamique [6.4] a relié le MV à la classe d'universalité du modèle XY. Récemment, une version discrétisée sur un réseau bidimensionnel a été étudiée, le modèle d'Ising actif (MIA) [6.5,6.6], pour lequel les particules actives sautent vers la gauche ou la droite avec un biais, et s'alignent localement comme des spins d'Ising. La transition vers un mouvement collectif est très similaire au MV, par contre ce modèle permet d'obtenir de nombreux résultats analytiques.

6.1 Mouvement collectif des spins de Potts

Nous avons considéré une généralisation pour le modèle de Potts à $q$ états. Les particules ont $q$ états internes correspondant à leurs directions de mouvement, et s'alignent localement comme des spins de Potts. Ce modèle est réaliste uniquement pour $q=4$ sur un réseau carré et $q=6$ sur un réseau triangulaire, $q=2$ revenant au MIA. Dans ce résumé, je n'évoquerai que les résultats pour $q=4$, reproductibles pour $q=6$. Considérons alors $N$ particules sur un réseau carré périodique doté de $L^2$ sites, pour une densité moyenne $\rho_0=N/L^2$. L'état de spin de la $k^{\rm e}$ particule sur le site $i$ est notée $\sigma_i^k$ avec une valeur entière comprise entre $0$ et $3$. Le nombre de particules dans l'état $\sigma$ sur le site $i$ est noté $n_i^\sigma$ tandis que la densité locale de ce site est $\rho_i = \sum_{\sigma=0}^3 n_i^\sigma$, sans aucune restriction sur sa valeur maximale. Sur le site $i$, la magnétisation est égale à $m_i^\sigma = (4n_i^\sigma -\rho_i)/3$ et l'hamiltonien est défini intensivement par $$H_i = - \frac{J}{2\rho_i} \sum_{k=1}^{\rho_i} \sum_{l \ne k} (q \delta_{\sigma_i^k,\sigma_i^l}-1), \tag{6.1}$$ avec une constante de couplage $J$ et le delta de Kronecker $\delta_{\alpha \beta}$. Pendant un intervalle de temps $\Delta t$, chaque particule située sur un site $i$ dans un état de spin $\sigma$ peut soit changer de spin vers l'état $\sigma'$ soit sauter vers un site voisin. La probabilité de ce changement de spin est $p_{\rm flip} = \exp [-4\beta J (n_i^\sigma - n_i^{\sigma'}-1) / \rho_i] \Delta t$, défini à partir de la différence d'énergie induite par l'hamiltonien $H_i$ pour une température $T=\beta^{-1}$. Et la probabilité de ce saut est $p_{\rm saut} = D(1+\epsilon) \Delta t$ dans la direction $\phi = \pi \sigma/2$ donnée par l'état de spin, et $p_{\rm saut} = D(1-\epsilon/3)\Delta t$ dans les trois autres directions. La probabilité totale de saut est $4D\Delta t$, indépendante du biais $\epsilon$. Des résultats numériques sont alors accessibles en utilisant l'algorithme de Metropolis [6.7].

Dans la limite hydrodynamique, à l'échelle mésoscopique ($L\gg1$), la densité de particules dans l'état $\sigma$ est définie par $\rho_\sigma({\bf x},t) = \langle n_i^\sigma \rangle$ et la densité totale de particules par $\rho({\bf x},t) = \langle \rho_i \rangle$, pour la coordonnée cartésienne ${\bf x}=(x,y)$ dont les positions entières correspondent aux nœuds du réseau. À partir de l'équation maîtresse du modèle microscopique, nous avons dérivé l'équation hydrodynamique $$\partial_t \rho_\sigma = D_\parallel \partial_\parallel^2 \rho_\sigma + D_\perp \partial_\perp^2 \rho_\sigma - v \partial_\parallel \rho_\sigma + \sum_{\sigma' \ne \sigma } \left[\frac{4\beta J}{\rho}(\rho_\sigma+\rho_{\sigma'}) -1 - \frac{r}{\rho} - \alpha \frac{(\rho_\sigma-\rho_{\sigma'})^2}{\rho^2}\right](\rho_\sigma-\rho_{\sigma'}), \tag{6.2}$$ avec une diffusivité longitudinale $D_\parallel = D(1+\epsilon/3)$ dans la direction ${\bf e_\parallel} = (\cos \phi, \sin \phi)$, une diffusivité transverse $D_\perp = D(1-\epsilon/3)$ dans la direction ${\bf e_\perp} = (-\sin \phi,\cos \phi)$ et une vitesse d'autopropulsion $v=4D\epsilon/3$. Les dérivées partielles sont définies par $\partial_\parallel = {\bf e_\parallel} \cdot \nabla$ et $\partial_\perp = {\bf e_\perp} \cdot \nabla$. Le dernier terme correspond aux changements de spins, dérivé en considérant que les magnétisations $m_i^\sigma$ sont distribuées identiquement avec une moyenne $\langle m_i^\sigma \rangle \ll \langle \rho_i \rangle$ et une variance proportionnelle à $\langle \rho_i \rangle$. Cette hypothèse cruciale a été vérifiée avec les simulations numériques. Les constantes $\alpha = 8 (\beta J)^2 (1-2\beta J/3)$ et $r$ dépendent uniquement de la température, et $r$ redimensionne les densités ce qui permet le choix $r=1$. Notons que la valeur $r=0$ correspond à l'approximation de champ moyen, sans prendre en compte les fluctuations, ne produisant que des états stationnaires homogènes. Toutes ces démarches sont similaires à celles utilisées pour développer l'équation hydrodynamique pour le MIA [6.5,6.6].

En redimensionnant le temps et la température par le choix $D=1$ et $J=1$, nous avons résolu numériquement l'Eq. (6.2) par la méthode des éléments finis avec FreeFem++ [6.8], dans un domaine carré de taille $L=100$ et en temps discrétisé. Nous nous sommes principalement intéressés à l'état stationnaire. Quand le biais est nul ($\epsilon=0$), une transition de phase a lieu pour une densité $\rho_0 = \rho_*$ à température fixée, causée par la brisure spontanée de la symétrie $Z_4$ des états de spin. Elle sépare une phase désordonnée à basse densité, pour laquelle tous les états sont identiquement peuplés, d'une phase ordonnée à haute densité, pour laquelle un état est surpeuplé par rapport aux trois autres identiquement peuplés. En augmentant la valeur du biais de $\epsilon=0$ à $\epsilon=3$ tout en gardant $\beta$ et $\rho_0=\rho_*$ fixes, nous observons une séparation de phase sous forme de bandes et d'allées, respectivement pour des petits et grands biais. Sur la Fig. 5a, une bande transverse au mouvement est tracée pour $\epsilon=1$ alors que sur la Fig. 5b, une allée longitudinale au mouvement est représentée pour $\epsilon=2.5$ pour une température inverse $\beta=0.75$ et une densité $\rho_0=1.33$. La phase ordonnée (liquide) se déplace sur la phase désordonnée (gaz) immobile. La réorientation d'une bande en une allée se produit pour un biais $\epsilon \simeq 2.0$ à cette température. Cette réorientation doit son origine à la diminution de la diffusion transverse $D_\perp$ et d'une augmentation de la vitesse $v$ quand le biais devient grand, stabilisant les allées. Il est à noter que cette réorientation est absente du MV et du MIA. À partir de ces résultats, nous avons tracé les diagrammes de phase, température-densité et biais-densité, reproduits sur les Fig. 55c et Fig. 5d. Comme pour le MV et le MIA, la transition vers le mouvement collectif est donc similaire à une transition de phase liquide-gaz sans région supercritique, avec des binodaux $\rho_g$ et $\rho_l$ correspondant aux densités des phases gazeuse et liquide dans la région de coexistence. Tous ces résultats sont concordants avec ceux rapportés par les simulations numériques.

Pour finir, nous avons dérivé les solutions stationnaires homogènes de l'Eq. (6.2) et nous avons étudié leur stabilité par une analyse linéaire, afin d'expliquer analytiquement la réorientation des bandes en allées. Pour une densité $\rho_0$ fixée, trois solutions homogènes différentes sont obtenues. Une solution pour la phase désordonnée avec une magnétisation $m_0=0$, et deux solutions pour la phase ordonnée avec une magnétisation $$\frac{m_0}{\rho_0}= \frac{\beta J}{\alpha} \left\{ 1 \pm \sqrt{1+\frac{(\beta-1-r/\rho_0)\alpha}{(\beta J)^2}} \right\} = \frac{\beta J}{\alpha} \pm \sqrt\frac{r(\rho_0-\rho_*)}{\alpha \rho_0\rho_*}, \tag{6.3}$$ définie pour une densité supérieure à $\rho_*$ et une température inférieure à $T_c = (1-\sqrt{22}/8)^{-1} \simeq 2.417$, cohérentes avec les simulations numériques. Au-dessus de cette température critique, seule la phase gazeuse existe. La magnétisation étant discontinue en $\rho_0=\rho_*$, la transition de phase pour $\epsilon=0$ est du premier ordre. Ces différentes valeurs de la magnétisation sont tracées sur la Fig. 5e. Nous avons alors analysé la stabilité de ces solutions en ajoutant une petite perturbation à celles-ci. L'évolution de cette perturbation est obtenue en linéarisant l'Eq. (6.2) et en calculant les valeurs propres de son opérateur d'évolution. Si une perturbation dans la direction longitudinale (resp. transverse) au mouvement s'amplifie alors une bande (resp. allée) se forme. Sur la Fig. 5f, nous avons tracé le diagramme de stabilité des bandes et allées correspondant pour $\beta=0.75$. Avec ces résultats, nous avons dérivé une expression analytique du biais pour lequel la réorientation a lieu dans la limite $T \rightarrow T_c$.

Nous avons donc décrit le mouvement collectif de spins de Potts possédant quatre et six états. La transition vers celui-ci est similaire à une transition de phase de type liquide-gaz sans région supercritique, comme pour le MIA (deux états) et le MV. De plus, nous avons observé une réorientation des bandes en allées, absente des modèles précédents, et expliqué son origine. Dans les prochains mois, nous allons étudier ce problème pour une autre interaction ferromagnétique, le modèle d'horloge allant vers le modèle XY lorsque le nombre d'état devient infini, pour comprendre les différences observées avec le MV. De plus, nous allons aborder ce modèle de Potts actif en ajoutant une restriction sur le nombre maximal de particules sur un site.

6.2 Mouvement collectif des spins de Potts en présence de restriction

6.3 Du modèle d'Ising au modèle de Vicsek

---

[P13] S. Chatterjee, M. Mangeat, R. Paul, and H. Rieger, Flocking and reorientation transition in the 4-state active Potts model, EPL 130, 66001 (2020).[P14] M. Mangeat, S. Chatterjee, R. Paul, and H. Rieger, Flocking with a q-fold discrete symmetry: Band-to-lane transition in the active Potts model, Phys. Rev. E 102, 042601 (2020).[P15] S. Chatterjee, M. Mangeat, and H. Rieger, Polar flocks with discretized directions: the active clock model approaching the Vicsek model, EPL 138, 41001 (2022).[P16] M. Karmakar, S. Chatterjee, M. Mangeat, H. Rieger, and R. Paul, Jamming and flocking in the restricted active Potts model, Phys. Rev. E 108, 014604 (2023).[6.1] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrodynamics of soft active matter, Rev. Mod. Phys. 85, 1143-1189 (2013).[6.2] T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and O. Shochet, Novel Type of Phase Transition in a System of Self-Driven Particles, Phys. Rev. Lett. 75, 1226 (1995).[6.3] A. P. Solon, H. Chaté, and J. Tailleur, From Phase to Microphase Separation in Flocking Models: The Essential Role of Nonequilibrium Fluctuations, Phys. Rev. Lett. 114, 068101 (2015).[6.4] J. Toner and Y. Tu, Long-Range Order in a Two-Dimensional Dynamical XY Model: How Birds Fly Together, Phys. Rev. Lett. 75, 4326 (1995).[6.5] A. P. Solon and J. Tailleur, Revisiting the Flocking Transition Using Active Spins, Phys. Rev. Lett. 111, 078101 (2013).[6.6] A. P. Solon and J. Tailleur, Pattern formation in flocking models: A hydrodynamic description, Phys. Rev. E 92, 042119 (2015).[6.7] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, Equation of state calculations by fast computing machines, J. Chem. Phys. 21, 1087–1092 (1953).[6.8] F. Hecht, New development in freefem++, J. Num. Math. 20, 251–266 (2013).

7. Modèle de Vicsek à deux espèces

---

[P17] S. Chatterjee, M. Mangeat, C.-U. Woo, H. Rieger, and J. D. Noh, Flocking of two unfriendly species: The two-species Vicsek model, Phys. Rev. E 107, 024607 (2023).

8. Sedimentation de particles actives

---

[P18] M. Mangeat, S. Chakraborty, A. Wysocki, and H. Rieger, Stationary particle currents in sedimenting active matter wetting a wall, Phys. Rev. E 109, 014616 (2024).