Keywords: tracking, superdiffusive transients, media, constraint, 2, 6 Extended crowding, gating functions

Introduction 2 Diffusion in crowded environments: from microscopic stochastic processes to modified diffusion models 2.1 The standard diffusion equation and its microscopic derivation . . . 2.2 Diffusion of tagged particles in a crowded medium . . . . . . . . . . 2.2

Chapter 1 Introduction

Diffusive transport is central in many areas of physics, chemistry, biology and soft matter [START_REF] Crank | The mathematics of diffusion[END_REF][START_REF] Dhont | An introduction to dynamics of colloids[END_REF][START_REF] Mcguffee | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF][START_REF] Caspi | Diffusion and directed motion in cellular transport[END_REF]. The regular functionality and development of cells, and in general of all biological organisms, could not prescind from the efficient delivery of a great number of molecules to their proper destination. New products are continuously synthetized and assembled inside the nucleus, and from there they have to reach different compartments. Conversely, degraded elements have to find their way back to the nucleus to be processed. Intracellular transport consists of two mechanisms: passive diffusion and active transport driven by molecular motors. Passive diffusion involves the movement of molecules through membranes (filtration) or within the cytosol, resulting from concentration gradients between two regions. Conversely, active transport forces particles to move against their concentration gradient. For this reason, it requires chemical energy. Living cells are indeed complex systems, whose maintenance involves a huge number of processes in which enzymes catalyze reactions, interact with each other and with cell products, and where proteins feedback to regulate the synthesis of other molecules and to trigger defensive measures in case of attack by hostile agents. In general, molecules diffuse in the hosting medium, chasing the target partner for a chemical reaction to take place. Thus, to the extent of providing an overall description of the cell metabolism, it is a matter of fundamental importance to characterize the steps through which the actors of the different reactions get in contact before binding. Moreover, dysfunctions in molecule transport can compromise seriously the integrity of biological tissues, as it has been recently pointed out by numerous studies which connected anomalous diffusion of proteins to neurodegenerative diseases such as Alzheimers [START_REF] Bressloff | Stochastic models of intracellular transport[END_REF]. Diffusion is a key process also for biomedical and industrial application and materials science. Catalysis, drug delivery, filtration, development of nanotechnologies are just some examples of activities where transport phenomena involving diffusion are central and call for a self-consistent description.

The classical approach to diffusion, based on Fick's law is strictly valid only for ideal, infinitely diluted solutions, where one is allowed to neglect the interaction among diffusing particles and between the molecules and the environment. Real media however, either biological structures or artificial technologies, are far from the diluted limit: cells for example contain a large number of membranes, compartments, and other smaller organelles that occupy up to 30 -40% of the available volume and that cannot overlap with each other. Crowding and confinement effects on diffusion-influenced phenomena pose hence fundamental yet unanswered questions. Concerning molecular mobility, for example, several computational and experimental indications exist of anomalous diffusion in the cell cytoplasm depending on the amount and type of crowding [START_REF] Balbo | The shape of protein crowders is a major determinant of protein diffusion[END_REF][START_REF] Weiss | Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells[END_REF], suggesting that living cells behave much like fractal or otherwise disordered systems [START_REF] Daniel Ben | Diffusion and Reactions in Fractals and Disordered Systems[END_REF][START_REF] Bouchaud | Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications[END_REF]. However, strong evidences also exist in favour of normal (Brownian) diffusion: in this scenario crowding and confinement result in (often nontrivial) modifications of the diffusion coefficient [START_REF] Dix | Crowding effects on diffusion in solutions and cells[END_REF][START_REF] Novak | Diffusion in cytoplasm: Effects of excluded volume due to internal membranes and cytoskeletal structures[END_REF][START_REF] Michael | Crowding and confinement effects on protein diffusion in vivo[END_REF][START_REF] Nicholson | Ion diffusion modified by tortuosity and volume fraction in the extracellular microenvironment of the rat cerebellum[END_REF]. Another related issue is that of diffusion-limited reactions [START_REF] Smoluchowski | Drei vortrage ubër diffusion brownsche molekular bewegung und koagulation von kolloidteichen[END_REF], which are ubiquitous in many domains of biology and chemistry, touching upon problems such as association, folding and stability of proteins [START_REF] Zhou | Protein folding and binding in confined spaces and in crowded solutions[END_REF][START_REF] Cheung | Molecular crowding enhances native state stability and refolding rates of globular proteins[END_REF] and bimolecular reactions in solution [START_REF] Dzubiella | Substrate concentration dependence of the diffusion-controlled steady-state rate constant[END_REF][START_REF] Piazza | Anticooperativity in diffusion-controlled reactions with pairs of anisotropic domains: A model for the antigen-antibody encounter[END_REF][START_REF] Fanelli | Diffusion in a crowded environment[END_REF][START_REF] Schmit | Lattice model of diffusion-limited bimolecular chemical reactions in confined environments[END_REF][START_REF] Tachiya | Concentration dependence of fluorescence quenching by ionic reactants[END_REF], including enzyme kinetics [START_REF] Agrawal | Effect of macromolecular crowding on the rate of diffusion-limited enzymatic reaction[END_REF], but also the dynamics of active agents [START_REF] Ricci | A quantitative approach to movement, displacement, and mobility of protozoa[END_REF]. Many theoretical models of intracellular diffusion have been derived in the past 30 years, and many aspects such as viscosity, crowding, confinement, and signaling mechanisms have been recently investigated in order to assess their impact on the transport of particles and the consequent reactions between components. Nevertheless, a full theoretical comprehension of transport in nonideal media remains an elusive task, Fick's law itself and the very notion of effective diffusion coefficient being questionable in a disordered medium [START_REF] Bouchaud | Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications[END_REF]. The aim of this thesis is to shed some light on the effect of crowding and complex geometries on the mobility of molecules, and to provide different approaches to include these aspects in the general investigation of diffusive processes and diffusion influenced reactions. In the first part of this work we will concentrate on the effect of crowding, describing at the microscopic level the interaction between diffusive agents, while the analysis of diffusion-regulated processes occurring in complex geometries will be the subject of the second part of the thesis.

The second chapter is devoted to characterize the influence of a densely populate medium on the transport of molecules, in different settings which account respectively for unbiased and biased diffusion in homogeneous and non homogeneous media. For all the studied problems we aim at deriving continuum models starting from microscopic stochastic processes with exclusion, which describe the dynamics of the single agents. We will show how the effect of crowding is lost when passing from the point-like agent based models to the continuum equations in the absence of bias and spatial anisotropy. Conversely, we will see that the excluded volume effects are always recognizable in the mean field models if some degree of asymmetry or inhomogeneity is present or if specific agents are tagged.

The third chapter deals with the analysis of diffusion-influenced reactions occurring between molecules which are diffusing in confined domains, accounting for the presence of obstacles and traps. Different degrees of chemical affinity between the reactants are envisioned. The general binding process is studied by solving the stationary diffusion equation using techniques of harmonic analysis, based on addition theorems for spherical harmonics [START_REF] Morse | Methods of Theoretical Physics, Part II[END_REF]. The model enables us to compute exactly the rate of a general pseudo first order reaction for an arbitrary configuration of bodies in a confined domain (the unbounded case is a simple generalization of this setting), and to analyze the diffusive interaction between particles that are reactive with the same species. Moreover, we provide closed formulas, obtained through simple approximations of the analytical solution and through perturbative expansions in powers of appropriate parameters, which allow one to characterize the depen-dence of the rate of the reaction on the relevant physico-chemical and geometrical parameters.

In chapters 4 and 5 we present two applications of our method. First, the analytical framework is employed within an industrial-chemical setting to analyze the catalytic optimization of a nanoreactor consisting of metal particles embedded in a thermosensitive microgel. Thereafter, we use the analytical method to evaluate the binding rate of an antibody-antigen reaction, and to identify the role of the protein's flexibility in the expression of the binding dynamics. Moreover we propose a simplified effective description which captures the essence of the investigated process.

In the last chapter, the time dependence is re-introduced to analyze gated diffusion-influenced reactions, where the binding process between a protein and a substrate is limited by structural fluctuations of the protein conformation which control the accessibility of the binding sites. The problem has been tackled through a suitable adaptation of Duhamel's theory to the diffusion equation with timedependent boundary conditions for an isolated receptor. The original problem of computing the rate of the gated reaction is thus reduced to the solution of a linear Volterra equation of the second kind, which is then solved with an appropriate algorithm. Finally, we present a brief analysis of the different regimes that can be classified depending on the function that regulates the accessibility of the site.

Chapter 2

Diffusion in crowded environments: from microscopic stochastic processes to modified diffusion models

In this chapter we will provide several models for molecular diffusion in crowded environments. As a common strategy for all the processes, we will base our analysis on the description of the motion of the individual agents, following the procedure presented in the next section.

The standard diffusion equation and its microscopic derivation

Diffusion is a fundamental process in nature that describes the spread of particles subject to random forces from regions of high density to regions of low density [START_REF] Crank | The mathematics of diffusion[END_REF]. The standard diffusion equation was originally derived by Adolf Fick in 1855. It is obtained by the combination of the first Fick's law:

J = -D∇ρ
which relates the diffusive flux J to the concentration ρ of particles, together with his second law, which is nothing but a continuity equation for the concentration:

∂ρ ∂t + ∇ • J = 0,
which states the conservation of the mass. The result, for a constant diffusion coefficient D, is the well-known diffusion equation:

∂ρ ∂t = D ∇ 2 ρ. (2.1)
However, while the mathematics of diffusive processes in dilute and simple media is fairly well developed and understood [START_REF] Crank | The mathematics of diffusion[END_REF], many interesting and relevant diffusive Diffusion in crowded environments processes take place in strongly non-ideal conditions. These include a wealth of different highly dense media, from non-ideal plasmas [START_REF] Kremp | Reaction and diffusion in dense nonideal plasmas[END_REF] to biological membranes [START_REF] Saffman | Brownian motion in biological membranes[END_REF], media with complex topological structures, including porous media [START_REF] Talmon | Molecular diffusion in porous media by pgse esr[END_REF][START_REF] Mitra | Diffusion propagator as a probe of the structure of porous media[END_REF][START_REF] Whitaker | Diffusion and dispersion in porous media[END_REF] and living cells [START_REF] Dix | Crowding effects on diffusion in solutions and cells[END_REF][START_REF] Phair | High mobility of proteins in the mammalian cell nucleus[END_REF] and strongly confining environments [START_REF] Condamin | Firstpassage times in complex scale-invariant media[END_REF][START_REF] Zhou | Protein folding and binding in confined spaces and in crowded solutions[END_REF][START_REF] Minton | Implications of macromolecular crowding for protein assembly[END_REF][START_REF] Mcguffee | Diffusion, crowding & protein stability in a dynamic molecular model of the bacterial cytoplasm[END_REF][START_REF] Ando | Crowding and hydrodynamic interactions likely dominate in vivo macromolecular motion[END_REF].

In many cases it is extremely complicated to understand how to modify the standard models of transport in order to account for non-ideal conditions. For several situations it may prove simpler or more effective to describe a complex transport process (or a simple one occurring in a complex medium) at the microscopic level. This means that, instead of describing the evolution of global quantities, such as the density of particles, we concentrate on the behavior of individual molecules which are seen as agents capable of moving according to specific rules. With this strategy we obtain different stochastic equations (corresponding to different deviations from the standard conditions) for the variation over time of the probability to find the system in a certain state. The macroscopic picture can then be recovered straightforwardly as the continuum limit of the agent-based model. As we will see in the following sections, the transition from the microscopic to the continuum description is a delicate procedure that often causes the loss of important information. It is thus important to investigate how the two levels of description interface with each other.

The standard diffusion equation (2.1) is the classic example of a macroscopic transport equation which can be obtained starting from an agent based model. Let us focus on the motion of the individual diffusing particles. If we denote with P i (n) the probability that an agent is at site i on some discrete manifold in one dimension at time n∆t, a simple unbiased random walk corresponds to the update rule

P i (n + 1) = 1 2 [P i-1 (n) + P i+1 (n)] (2.2)
as it is assumed that at each time step the walker can either jump to its right or to its left with equal probability. Thus introducing the time step ∆t and the spacing ∆x we have:

P i (n + 1) -P i (n) ∆t = ∆x 2 ∆t 1 2∆x 2 [P i-1 (n) + P i+1 (n) -2P i (n)], (2.3)
where we can recognize the discrete laplacian on the right-hand side.

Letting the lattice spacing ∆x and the time step ∆t go to zero, such that lim ∆x,∆t→0

(∆x 2 /2∆t) = D, one obtains the diffusion equation (2.1) in the continuum limit.

In performing the continuum limit one is tacitly assuming that many walkers are performing as many uncorrelated random walks and that a probability of being at x at time t can be defined by averaging over such uncorrelated trajectories. This requires the walkers to be transparent to each other. It is interesting to ask the following question. If some exclusion rule is enforced, such that a walker can only jump on an empty site, how will the macroscopic equation be modified? And what kind of process will they describe?

Starting from microscopic models we want to understand to which extent relaxing the hypotheses of diluted medium, homogeneous space and absence of external potentials leads to a modification of the standard diffusion equation.

Diffusion of tagged particles in a crowded medium 2.2.1

Effects of crowding on the diffusion properties: normal or anomalous diffusion?

The study of molecular diffusion under crowded conditions represents a particularly crucial topic for its applications to cellular biology. At high density, particles diffusion is impeded and excluded-volume effects may no longer be ignored. The hallmark of diffusive transport is the linear growth in time of the mean square displacement (MSD) of the spreading particles: ∆R 2 ∝ t. This is a simple conclusion that follows directly from the law of conservation of matter (in the form of a continuity equation), when a simple constitutive equation is assumed, stating that the particle current is proportional to the concentration gradient. The latter law, known as the (first) Fick's law, can be regarded as a simple linear-response prescription, thus only appropriate to describe the relaxation of small density fluctuations.

Despite the fact that Fickean diffusion is generally appropriate to describe the spontaneous spatial rearrangement of particles in suspension, deviations are expected to occur in various situations of interest, e.g. if fixed obstacles are present (confinement) [START_REF] Sekhar Burada | Diffusion in confined geometries[END_REF] or when different, and thus distinguishable species compete for the available space at high concentration, a scenario often referred to in cellular biology as macromolecular crowding [START_REF] Foffi | Macromolecular crowding: chemistry and physics meet biology[END_REF][START_REF] Rivas | Life in a crowded world[END_REF][START_REF] Dix | Crowding effects on diffusion in solutions and cells[END_REF][START_REF] John | Macromolecular crowding: obvious but underappreciated[END_REF][START_REF] John | Macromolecular crowding: an important but neglected aspect of the intracellular environment[END_REF][START_REF] Michael | Crowding and confinement effects on protein diffusion in vivo[END_REF][START_REF] Allen P Minton | Influence of macromolecular crowding upon the stability and state of association of proteins: Predictions and observations[END_REF][START_REF] Schnell | Reaction kinetics in intracellular environments with macromolecular crowding: simulations and rate laws[END_REF][START_REF] Weiss | Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells[END_REF][START_REF] Zhou | Protein folding and binding in confined spaces and in crowded solutions[END_REF][START_REF] Zhou | Macromolecular crowding and confinement: biochemical, biophysical, and potential physiological consequences[END_REF][START_REF] Szymański | Diffusion and Viscosity in a Crowded Environment: from Nano-to Macroscale[END_REF][START_REF] Phair | High mobility of proteins in the mammalian cell nucleus[END_REF][START_REF] Piazza | The unbearable heaviness of colloids: facts, surprises, and puzzles in sedimentation[END_REF]. Despite the importance of crowding and confinement effects in diffusion-related mechanisms in chemistry and biology, there is no consensus on the mechanisms through which crowding and confinement fine-tune deviations from the classical Fickean picture. This lively debate is reflected by conflicting experimental reports in the literature concerning the role of complex environmental factors in the mobility of biomolecules in the cytoplasm and extra-cellular matrix. Some authors maintain that crowding merely slows down transport by reducing in a complex fashion the diffusion coefficient but does not alter the MSD exponent [START_REF] Dix | Crowding effects on diffusion in solutions and cells[END_REF][START_REF] Novak | Diffusion in cytoplasm: Effects of excluded volume due to internal membranes and cytoskeletal structures[END_REF][START_REF] Dauty | Molecular crowding reduces to a similar extent the diffusion of small solutes and macromolecules: measurement by fluorescence correlation spectroscopy[END_REF][START_REF] Szymański | Diffusion and Viscosity in a Crowded Environment: from Nano-to Macroscale[END_REF], while others [START_REF] Golding | Physical Nature of Bacterial Cytoplasm[END_REF][START_REF] Pastor | Diffusion of alpha-chymotrypsin in solutioncrowded media. A fluorescence recovery after photobleaching study[END_REF][START_REF] Weiss | Anomalous subdiffusion is a measure for cytoplasmic crowding in living cells[END_REF][START_REF] Daniel | Anomalous Diffusion of Proteins Due to Molecular Crowding[END_REF] contend the identification of crowding with anomalous (typically sub-diffusive) transport [START_REF] Bouchaud | Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications[END_REF][START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF][START_REF] Zaslavsky | Chaos, fractional kinetics, and anomalous transport[END_REF][START_REF] Denis S Grebenkov | Hydrodynamic and subdiffusive motion of tracers in a viscoelastic medium[END_REF], a feature observed in lateral diffusion in cellular membranes [START_REF] Javanainen | Anomalous and normal diffusion of proteins and lipids in crowded lipid membranes[END_REF][START_REF] Gerald R Kneller | Communication: Consistent picture of lateral subdiffusion in lipid bilayers: Molecular dynamics simulation and exact results[END_REF][START_REF] Feder | Constrained diffusion or immobile fraction on cell surfaces: a new interpretation[END_REF]. In this case one would have ∆R 2 ∝ t α with α < 1 (sub-diffusion) or ∆R 2 ∝ t α with α > 1 (super-diffusion). It is worthwhile to underline that reports of anomalous transport connected to crowding are not limited to sub-diffusion. For example, Upadhaya and collaborators [START_REF] Upadhyaya | Anomalous diffusion and non-Gaussian velocity distribution of Hydra cells in cellular aggregates[END_REF] have recorded super-diffusive behavior in the motion of endodermal Hydra cells, which they traced back to long-range correlations within the scrutinized sample, while Stauffer and collaborators [START_REF] Stauffer | Superdiffusion in a Model for Diffusion in a Molecularly Crowded Environment[END_REF] proposed a minimalistic model of random barriers in a percolation network as a tool to mimic diffusion in a crowded environment.

To add an important piece of information to the debate, it is interesting to remark that most often claims of anomalous diffusion in three-dimensional crowded environments in vitro and in vivo rely on fluorescence recovery after photobleaching (FRAP) data that are analyzed through ad hoc modifications [START_REF] Pastor | Diffusion of alpha-chymotrypsin in solutioncrowded media. A fluorescence recovery after photobleaching study[END_REF] of standard theories of fluorescence photobleaching recovery [START_REF] Soumpasis | Theoretical analysis of fluorescence photobleaching recovery experiments[END_REF][START_REF] Axelrod | Mobility measurement by analysis of fluorescence photobleaching recovery kinetics[END_REF]. It is interesting to remark that to our knowledge no first-principle derivations of fluorescence recovery curves in the anomalous diffusion regime have yet been reported, analogous to the long-known standard derivations performed in the context of normal diffusion [START_REF] Soumpasis | Theoretical analysis of fluorescence photobleaching recovery experiments[END_REF][START_REF] Axelrod | Mobility measurement by analysis of fluorescence photobleaching recovery kinetics[END_REF].

i i 1 i+1

Figure 2.1: The simple exclusion rule implemented in the model of tagged particle diffusion. In this configuration, the tagged particle (filled circle) sitting at site i can only jump towards the (empty) neighboring site i -1. The jump towards site i + 1 is instead impeded, as the target site is occupied by a crowder (empty circles). Crowders can also diffuse towards neighboring sites.

As it is often the case, the truth probably reflects an intermediate picture. Possibly, complex (even multiple) crossovers are to be expected between anomalous and normal diffusion [START_REF] Vilaseca | New insights into diffusion in 3D crowded media by Monte Carlo simulations: effect of size, mobility and spatial distribution of obstacles[END_REF], or, alternatively, one needs to consider complex spaceand geometry-dependent diffusion coefficients [START_REF] Zador | Microfiberoptic fluorescence photobleaching reveals size-dependent macromolecule diffusion in extracellular space deep in brain[END_REF][START_REF] Nicholson | Ion diffusion modified by tortuosity and volume fraction in the extracellular microenvironment of the rat cerebellum[END_REF], as modeled e.g. by Fick-Jacobs [START_REF] Jacobs | Diffusion Processes[END_REF] and related theories [START_REF] Martens | Transport of Brownian particles in confined geometries -Steps beyond the Fick-Jacobs approach[END_REF][START_REF] Kalinay | When is the next extending of Fick-Jacobs equation necessary[END_REF]. However, as it appears clear from the above recollection, the need for further, systematic investigation of transport in crowded and confining media is evident.

We will investigate hereafter the effect of crowding, by enforcing, in the microscopic model of a random walk, the exclusion rule which prevents the agent to jump to an occupied site. We will consider for the moment an homogeneous space, in the absence of external forces.

Microscopic processes implementing exclusion rules go under the name of simple exclusion processes (SEP). In general, SEPs are space-discrete, agent-based stochastic processes modeling some kind of transport according to specific rules and bound to the constraint that no two agents can ever occupy the same site. SEPs play a central role in non-equilibrium statistical physics [START_REF] V Privman | Nonequilibrium Statistical Mechanics in One Dimension[END_REF][START_REF] Liggett | Stochastic Interacting Systems: Contact, Voter and Exclusion Processes[END_REF]. While the first theoretical ideas underlying such processes can be traced back to Boltzmann's works [START_REF] Boltzmann | Vorlesungen über Gastheorie (Lectures on Gas Theory), volume I and II[END_REF], SEPs were introduced and widely studied in the 70s as simplified models of one-dimensional transport for phenomena like hopping conductivity [START_REF] Peter | Theory of one-dimensional hopping conductivity and diffusion[END_REF] and kinetics of biopolymerization [START_REF] Bressloff | Stochastic models of intracellular transport[END_REF]. Along the same lines, the asymmetric exclusion process (ASEP), originally introduced by Spitzer [START_REF] Spitzer | Interaction of Markov processes[END_REF], has become a paradigm in non-equilibrium statistical physics [START_REF] Derrida | Exact solution of a 1d asymmetric exclusion model using a matrix formulation[END_REF][START_REF] Schütz | Phase transitions in an exactly soluble onedimensional exclusion process[END_REF][START_REF] Derrida | An exactly soluble non-equilibrium system: The asymmetric simple exclusion process[END_REF] and has now found many applications, such as the study of molecular motors [START_REF] Golubeva | Efficiency at maximum power of interacting molecular machines[END_REF], transport in the cytoskeleton [START_REF] Parmeggiani | Totally asymmetric simple exclusion process with langmuir kinetics[END_REF] through nano-channels [START_REF] Zilman | Crowding effects in non-equilibrium transport through nano-channels[END_REF] and depolymerization of microtubules [START_REF] Reese | Crowding of molecular motors determines microtubule depolymerization[END_REF].

The idea is to move from a space-discrete simple exclusion process specifying the competition for space at the microscopic level: the probability of an agent jumping from its current position in i on a neighboring site (i ± 1) is conditioned to the probability (1 -P i±1 ) of finding the site empty. If we enforce the excluded volume effect in the master equation for a single species diffusing in an isotropic space we obtain again the standard diffusion equation in the continuum limit. The excluded volume constraint is lost in the continuum limit: that happens because the agents are not distinguishable. By contrast if we make the agent distinguishable, either accounting for different species or (as we will show hereafter) tagging individual particles belonging to the same population, we obtain a modified diffusive behavior. More specifically cross-diffusive terms linking multiple diffusive components appear, modifying the standard Laplacian term [START_REF] Fanelli | Diffusion of two molecular species in a crowded environment: theory and experiments[END_REF]. If one takes into account two populations, in [START_REF] Fanelli | Diffusion of two molecular species in a crowded environment: theory and experiments[END_REF] the following macroscopic equations are obtained for the evolution of their densities φ and ψ:

           ∂φ ∂t = D φ ∇ 2 φ + φ∇ 2 ψ -ψ∇ 2 φ ∂ψ ∂t = D ψ ∇ 2 ψ + ψ∇ 2 φ -φ∇ 2 ψ (2.4)
where the cross terms φ∇ 2 ψ -ψ∇ 2 φ and ψ∇ 2 φ -φ∇ 2 ψ derive from imposing the excluded volume constraints in the microscopic equations.

Another way to keep track of the excluded volume effect in the macroscopic dynamic is by making a few molecules recognizable inside a single population of diffusive agents and monitoring the evolution of the system. In the following sections we will consider the diffusive process of an ensemble of particles, the tagged species, immersed in a densely populated background of co-evolving agents, hereafter the crowders.

This is a quite general scenario, which can be invoked to describe different experimental conditions, typical of fluorescence-based single-molecule tracking experiments. Hereafter we will present briefly the main experimental techniques typically employed to track single molecules.

Experimental particle tracking

The direct observation of the motion of single particles in vivo is a powerful method for describing the interaction between a diffusing molecule and the surrounding environment. Numerous optical techniques have been developed to provide a characterization of the diffusion of particles into living cells, in terms of their velocity and diffusion coefficients. Single Particle Tracking (SPT) consists in observing single molecules, which are made distinguishable with different probes (such as fluorophores) that can be excited by laser light at specific intensities. The individual trajectories are retraced through images acquired at different times, in which the coordinates of the single molecules are recorded, and subsequently assembled to reproduce the motion. The main limitation of SPT is that it is not suitable for measuring diffusion in three dimensions (because of the high velocity of three dimensional diffusion and the fact that the images obtained are projected on a surface), thus this technique is used mainly for membrane diffusion.

For diffusion in the cytosol, 3D-FRAP (Fluorescence Recovery After Photobleaching) has proven to be a very useful technique. In FRAP a selected volume is at first visualized at low intensity light, and the fluorescence (arising because the monitored diffusing particles are fluorescently labeled) is recorded. Then the region Diffusion in crowded environments is invested with higher intensity light, causing the fluorophores to bleach and thus reducing the fluorescence. After the fluorescence in the area has reached the minimum value, it starts increasing again because of the unbleached molecules that are diffusing back inside the region. The analysis of the fluorescence intensity curve allows to obtain information on the diffusive parameters. The diffusion constant D can be calculated as

D = R 2 4t D
where R is the radius of the beam (the profile is assumed to be gaussian) and t D is the characteristic time of diffusion.

Derivation of the equations for the tagged particles and for the bulk species

We will introduce now a strategy to describe the diffusion of recognizable agents among a dense population. The tagged particles are assumed to be sufficiently diluted. Following a microscopic approach, we derive a system of partial differential equations for the mean-field densities of both the tagged particles and the crowders. The model is formulated at the microscopic level as a stochastic process with simple exclusion interference. In the thermodynamic limit, the excluded-volume constraints will result in nonlinear coupling terms between the two concentrations. To simplify the discussion, let us consider a one-dimensional problem. As shown in the following, the derivation can be readily extended to higher dimensions without altering the ensuing physical picture. Let us consider a one-dimensional lattice of spacing a.

Each site can be occupied by either a crowder or a tagged particle. We denote with the binary variables m i (k) and n i (k) the occupancy of site i at time t = k∆t for the tagged and crowding particles, respectively. Hence m i (k), n i (k) can be either zero or one depending on whether site i is occupied or not by the respective particle.

The stochastic process that governs jumps of the tagged particles can be cast in the following form

m i (k + 1) -m i (k) = (z + i-1 m i-1 (k) + z - i+1 m i+1 (k))[1 -m i (k)][1 -n i (k)] -z + i m i (k)[1 -m i+1 (k)][1 -n i+1 (k)] -z - i m i (k)[1 -m i-1 (k)][1 -n i-1 (k)] (2.5)
Eq. (2.5), and its analogue for species n i (k), can be regarded as the update rule for a simple Monte Carlo process. If the target site is occupied by either a crowder or a tagged particle, the move cannot occur. The quantities z ± i are variables that take the value 0 or 1 depending on a random number ξ i uniformly distributed between 0 and 1. By considering homogeneous jump probabilities, q ± j = q for j = i, i ± 1, one can formally write

z + i-1 = θ(ξ i ) -θ(ξ i -q) (2.6) z - i+1 = θ(ξ i -q) -θ(ξ i -2q) z + i = θ(ξ i -2q) -θ(ξ i -3q) z - i = θ(ξ i -3q) -θ(ξ i -4q)
where θ(x) is the Heaviside step function and we are assuming q ≤ 1/4. Eqs. (2.56) entail z ± j = q, where . . . denotes an average over many values of ξ i for a fixed configuration {n i , m i }. The above process is entirely determined by the jump probabilities q, which we here assume constant and homogeneous.

A (discrete-time) master equation can be obtained by averaging over many Monte Carlo realizations performed according to the rule (2.5) and starting from the same initial condition (we denote this average by . . . ). Introducing the one-body occupancy probabilities ρ i (k) = m i (k) and φ i (k) = n i (k) and assuming a mean-field factorization for the two-body and three-body correlations, one finds Diffusion in crowded environments

equations ρ i (k + 1) -ρ i (k) = q (ρ i-1 (k) + ρ i+1 (k)) [1 -ρ i (k)] [1 -φ i (k)] -q ρ i (k)[2 -(ρ i-1 (k) + ρ i+1 (k)) -(φ i-1 (k) + φ i+1 (k)) + φ i+1 (k)ρ i+1 (k) + φ i-1 (k)ρ i-1 (k)] φ i (k + 1) -φ i (k) = w (φ i-1 (k) + φ i+1 (k)) [1 -φ i (k)] [1 -ρ i (k)] -w φ i (k)[2 -(φ i-1 (k) + φ i+1 (k)) -(ρ i-1 (k) + ρ i+1 (k)) + φ i+1 (k)ρ i+1 (k) + φ i-1 (k)ρ i-1 (k)],
where w denotes the jump probability associated with crowders' motion. To proceed in the analysis, we assume that the concentration of tagged particles is small, ρ i ≪ 1. We therefore approximate the previous equations as

ρ i (k + 1) -ρ i (k) = q (ρ i-1 (k) + ρ i+1 (k)) [1 -φ i (k)] -qρ i (k)[2 -(φ i-1 (k) + φ i+1 (k))] φ i (k + 1) -φ i (k) = w (φ i-1 (k) + φ i+1 (k) -2φ i (k)) .
(2.7)

Eq. (2.7) tells us that the microscopic exclusion constraint is lost in the equation for the crowders occupancy probability φ i . The tagged particles are in fact highly diluted and thus they interfere negligibly with the diffusive motion of the crowders.

Let us now move to the continuum. This can be formally achieved by letting

ρ(x, t) = lim a,∆t→0 ρ i (k) and φ(x, t) = lim a,∆t→0 φ i (k),
where we recall that a is the lattice spacing and ∆t is the time step. In addition we must require lim a,∆t→0 qa 2 /∆t = D ρ and lim a,∆t→0

wa 2 /∆t = D φ ,
where D ρ and D φ denote the diffusion coefficients of the tagged particles and the crowders, respectively. Making use of the above definitions, one readily obtains the continuum limit of Eqs. (2.7)

∂φ ∂t = D φ ∂ 2 φ ∂x 2 ∂ρ ∂t = D ρ ∂ ∂x ∂ ∂x [(1 -φ)ρ] + 2ρ ∂φ ∂x . (2.8)
The mean-field density of crowders φ evolves in time following a standard diffusion equation: the undisturbed dynamic is a consequence of the hypothesis of dilution of the tagged particles, which allowed us to neglect their feedback on the crowders. On the contrary, the density ρ obeys a Smoluchowski equation where the density of crowders plays the role of an external potential, which embodies the excludedvolume rules imposed at the microscopic level. In the next section we give an alternative derivation of eqs.(2.8), following a perturbative calculation inspired by the Van Kampen system size expansion [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF]. We note that the equation for the evolution of ρ has also been derived in Ref. [START_REF] Seki | Diffusion-mediated geminate reactions under excluded volume interactions[END_REF] for a constant non-homogeneous background field φ(x).

A derivation analogous to the one described above can be repeated in higher dimensions (see the next section), leading to a straightforward generalization of Eqs. (2.8)

∂φ(x, t) ∂t = D φ ∇ 2 φ(x, t) ∂ρ(x, t) ∂t + ∇ • J(x, t) = 0 (2.9)
where J = -D ρ {∇[(1φ)ρ] + 2ρ∇φ} is the total (osmotic plus force) current for the tagged species. In section 2.2.6 the above equations are integrated numerically, both in one and two dimensions. As we shall demonstrate, the effective force term leads to the emergence of sub-diffusive or super-diffusive transients in the dynamics of the tagged species, depending on the chosen initial conditions.

Derivation of the mean-field equations in higher dimensions and alternative approaches

As stated before, the equations governing the evolution of the particles, which have been obtained for a 1D setting, can be extended to higher dimensions. In the following we will provide a detailed derivation for the mean field equations in 2D.

We will first generalize straightforwardly to a higher dimensional space the procedure described in the previous section and then turn to consider two alternative approaches: the first is inspired to the work of Landman and collaborators (see e.g. [START_REF] Simpson | Pathlines in exclusion processes[END_REF]), while the second makes use of the Van Kampen expansion [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF].

To progress in the analysis we assume each site of the two-dimensional lattice to be labeled with two indices (i, j). In the 2D case, the selected particle can jump to four nearest neighboring sites. This is a slight modification of the one-dimensional geometry where each site counts only two adjacent neighbors. The binary variables at time k are hence here labeled with m i,j (k) and n i,j (k).

The stochastic process reads:

m i,j (k + 1) -m i,j (k) = z + i-1 m i-1,j (k)[1 -m i,j (k)][1 -n i,j (k)] + z - i+1 m i+1,j (k)[1 -m i,j (k)][1 -n i,j (k)] + z + j-1 m i,j-1 (k)[1 -m i,j (k)][1 -n i,j (k)] + z - j+1 m i,j+1 (k)[1 -m i,j (k)][1 -n i,j (k)] -z + i m i,j (k)[1 -m i+1,j (k)][1 -n i+1,j (k)] -z - i m i,j (k)[1 -m i-1,j (k)][1 -n i-1,j (k)] -z + i m i,j (k)[1 -m i,j+1 (k)][1 -n i,j+1 (k)] -z - i m i,j (k)[1 -m i,j-1 (k)][1 -n i,j-1 (k)],
with the stochastic variables z ± defined in analogy to the one dimensional case. The equation governing the evolution of n i,j (•) can be equivalently modified. After introducing the one-body occupancy probabilities

ρ i,j (k) = m i,j (k) (2.10) φ i,j (k) = n i,j (k) (2.11)
and assuming a mean-field factorization for the two-body and three-body correlations (the validity of these assumptions will be checked in Sec. 2.2.5 ), one eventually Diffusion in crowded environments ends up with

ρ i,j (k + 1) -ρ i,j (k) = q (ρ i-1,j (k) + ρ i+1,j (k) + ρ i,j-1 (k) + ρ i,j+1 (k)) [1 -ρ i,j (k)] [1 -φ i,j (k)] -q ρ i,j (k)[4 -(ρ i-1,j (k) + ρ i+1,j (k)) -(ρ i,j-1 (k) + ρ i,j+1 (k)) -(φ i-1,j (k) + φ i+1,j (k)) -(φ i,j-1 (k) + φ i,j+1 (k + 1)) + φ i+1,j (k)ρ i+1,j (k) + φ i-1,j (k)ρ i-1,j (k) + φ i,j+1 (k)ρ i,j+1 (k) + φ i,j-1 (k)ρ i,j-1 (k)] φ i,j (k + 1) -φ i,j (k) = w (φ i-1,j (k) + φ i+1,j (k) + φ i,j-1 (k) + φ i,j+1 (k)) [1 -φ i,j (k)] [1 -ρ i,j (k)] -w φ i,j (k)[4 -(φ i-1,j (k) + φ i+1,j (k)) -(φ i,j-1 (k) + φ i,j+1 (k)) -(ρ i-1,j (k) + ρ i+1,j (k)) -(ρ i,j-1 (k) + ρ i,j+1 (k)) + φ i+1,j (k)ρ i+1,j (k) + φ i-1,j (k)ρ i-1,j (k) + φ i,j+1 (k)ρ i,j+1 (k) + φ i,j-1 (k)ρ i,j-1 (k)].
If one assumes the concentration of the tagged particles to be small, namely ρ i,j ≪ 1, then the following approximated relations are found:

ρ i,j (k + 1) -ρ i,j (k) = q(ρ i-1,j (k) + ρ i+1,j (k) + ρ i,j-1 (k) + ρ i,j+1 (k)) [1 -φ i,j (k)] -qρ i,j (k)[4 -(φ i-1,j (k) + φ i+1,j (k) + φ i,j-1 (k) + φ i,j+1 (k))] φ i,j (k + 1) -φ i,j (k) = w(φ i-1,j (k) + φ i+1,j (k) -2φ i,j (k) + φ i,j-1 (k) + φ i,j+1 (k) -2φ i,j (k))
After introducing the continuous variables

ρ(x, y, t) = lim a,∆t→0 ρ i,j (k), φ(x, y, t) = lim a,∆t→0 φ i,j (k), (2.12)
where a and∆t respectively stand for the linear size of the lattice site and the characteristic time step of the microscopic dynamics and after defining the diffusion coefficients as it has been done for the one dimensional case, we get the sought-for generalized model:

∂φ ∂t = D φ ∇ 2 φ (2.13) ∂ρ ∂t = ∇ 2 [D ρ (1 -φ)ρ] + 2D ρ ∇ • (ρ∇φ) .
As anticipated, we shall now turn to discussing an alternative derivation of the above equations, obtained in line with the procedure used in [START_REF] Landman | Myopic random walkers and exclusion processes: Single and multispecies[END_REF]. To this end, we define the variable

γ i,j (k) =      1, if site (i, j
) is occupied by an agent of species m at time step k 2, if site (i, j) is occupied by an agent of species n at time step k 0 if site (i, j) is empty at the time step k

(2.14)
We write then the master equation for the motion of the tagged particles as:

P 1 (γ i,j (k + 1) = 1) -P 1 (γ i,j (k) = 1) = -α P 2 γ i,j (k) = 1, γ i+1,j (k) = 0 + P 2 γ i,j (k) = 1, γ i-1,j (k) = 0 + P 2 γ i,j (k) = 1, γ i,j+1 (k) = 0 + P 2 γ i,j (k) = 1, γ i,j-1 (k) = 0 + α P 2 γ i,j (k) = 0, γ i+1,j (k) = 1 + P 2 γ i,j (k) = 0, γ i-1,j (k) = 1 + P 2 γ i,j (k) = 0, γ i,j+1 (k) = 1 + P 2 γ i,j (k) = 0, γ i,j-1 (k) = 1 (2.15)
where P 2 stands for a joint probability, while P 1 is the probability of a single event and α is the rate of success of the selected jump. The master equation for the population n is similar, the variable γ assuming the value 2 instead of 1 and β labeling the associated jump rate (in principle the two population can have different jump probabilities). In the mean field limit, we factorize the joint probabilities P 2 in the master equations as:

P 2 γ i,j (k) = 1, γ i,j+1 (k) = 0 = P 1 γ i,j (k) = 1 P 1 γ i,j+1 (k) = 0 .
To perform the continuum limit we use the following Taylor expansion for the three different probability functions P 1 :

P(γ i±1,j (k) = 1) = ρ ± a ∂ρ ∂x + 1 2 a 2 ∂ 2 ρ ∂x 2 + o(a 2 ) P(γ i,j±1 (k) = 1) = ρ ± a ∂ρ ∂y + 1 2 a 2 ∂ 2 ρ ∂y 2 + o(a 2 ) P(γ i±1,j (k) = 2) = φ ± a ∂φ ∂x + 1 2 a 2 ∂ 2 φ ∂x 2 + o(a 2 ) P(γ i,j±1 (k) = 2) = φ ± a ∂φ ∂y + 1 2 a 2 ∂ 2 φ ∂y 2 + o(a 2 ) P(γ i±1,j (k) = 0) = µ ± a ∂µ ∂x + 1 2 a 2 ∂ 2 µ ∂x 2 + o(a 2 ) P(γ i,j±1 (k) = 0) = µ ± a ∂µ ∂y + 1 2 a 2 ∂ 2 µ ∂y 2 + o(a 2 ).
(2.16)

where, as stated before, a represents the linear size of each site. Making use of the relation µ = 1 -ρφ and defining the usual limits: lim a,∆t→0

αa 2 ∆t = D ρ lim a,∆t→0 βa 2 ∆t = D φ ,
we obtain the following equation for the tagged species:

∂ρ ∂t = D ρ ((1 -φ)∇ 2 ρ + ρ∇ 2 φ).

Diffusion in crowded environments

Under the hypotesis of low concentration of the tagged agents, we get the coupled equation for the evolution of the bulk density:

∂φ ∂t = D φ ∇ 2 φ.
We want to introduce now another alternative approach that leads to the continuum model (2.8), through a coarse-grained description of the scrutinized problem. The procedure is here explained in one dimension, but it readily generalizes to the relevant d = 3 setting. We consider the physical space to be partitioned in Ω patches, also called urns. Each patch has a maximum carrying capacity (it can be filled with N particles at most). Labeling m i the number of tagged particles contained in urn i, and with n i the corresponding number of crowders, one can write:

n i + m i + v i = N ∀i
where v i stands for the number of vacancies, the empty cases in patch i that can be eventually filled by incoming particles. The excluded-volume effect is here enforced by requiring that the particles can move only into the nearest-neighbor patches that exhibit vacancies, as exemplified by the following chemical reactions:

M i + V j δ zΩ -→ M j + V i N i + V j δ zΩ -→ N j + V i .
(2.17)

Here z is the number of nearest-neighbor patches and M i , N i , V i are respectively a particle of type M (the tagged particles), of type N (the crowders) or a vacancy belonging to the i-patch. This is a stochastic process which is governed, under the Markov hypothesis, by a Master equation for the probability P (n, m, t) of finding the system in a given state which is specified at time t by the 2Ω dimensional vector (n, m) = (n 1 , ..., n Ω , m 1 ..., m Ω ). The Master equation reads:

∂P (n, m, t) ∂t = n =n ′ ,m =m ′ [T (n, m|n ′ , m)P (n ′ , m) + T (n, m|n, m ′ )P (n, m ′ ) (2.18) -T (n, m ′ |n, m)P (n, m) -T (n ′ , m|n, m)P (n, m)],
where T (a|b) is the rate of transition from a state b to a compatible configuration a. The allowed transitions are those that take place between neighboring patches as dictated by the chemical reactions (2.17). For example, the transition probability associated with the second of equations (2.17) reads

T (n i -1, n j + 1|n i , n j ) = δ zΩ n i N v j N = δ zΩ n i N 1 - n j N - m j N . (2.19)
The transition rates introduce to the equation the explicit dependence on the system size, that is represented by the number N of molecules per patch. To proceed in the analysis, we adopt the Van Kampen system size expansion [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF], which enables us to separate the site-dependent mean concentration φ i (t) from the corresponding fluctuations ξ i in the expression of the discrete number density of species N . The role of fluctuations becomes less important as the number of the agents is increased, as made clear by the following Van Kampen ansatz:

n i N (t) = φ i (t) + ξ i √ N .
(2.20)

In the following we will assume the presence of a single tagged particle. However, the analysis extends straightforwardly to the case where more diluted particles are dispersed in the background of crowders. Since the tagged particle must be located in only one of the patches, it is convenient to look at the evolution of the function:

P k (n, t) = P (n, 0, 0, ..., 0 k-1
, 1, 0, ...., 0, 0

Ω-k , t)
inside the master equation (2.18). P k (n, t) is the probability that the target particle is found in the k-patch, for a particular configuration n of species N . With these notation the Master equation can be written in the following compact form:

∂P k (n, t) ∂t = Ω i=1 j∈i-1,i+1 (ǫ - j ǫ + i -1)T (n i -1, n j + 1|n i , n j )P k (n i , n j , t) + Ω i=1    - j∈i-1,i+1 δ zΩ 1 N 1 - n j N P k + j∈i-1,i+1 δ zΩ 1 N 1 - n k N P j    (2.21)
that involves the shift operators:

ǫ ± i f (...., n i , .....) = f (...., n i ± 1, .....).
Using the Van Kampen hypothesis [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry[END_REF], one can expand the transition rates in power of 1/ √ N . As a result, equation (2.19) takes the form

T (n i -1, n j + 1|n i , n j ) = δ zΩ (φ i (1 -φ j )) + 1 √ N ξ i (1 -φ j ) -ξ j φ i ] + 1 N [-ξ i ξ j -m j φ i ] + 1 N 3 2 -m j ξ i .
Moreover, the shift operators can be expressed as differential operators:

(ǫ - j ǫ + i -1) = 1 √ N ∂ ∂ξ i - ∂ ∂ξ j + 1 2N ∂ ∂ξ i - ∂ ∂ξ j 2 + O( 1 N 3 2
) so that

(ǫ - j ǫ + i -1)T (n i -1, n j + 1|n i , n j )P (n i , n j , t) = 1 √ N ∂ ∂ξ i - ∂ ∂ξ j δ zΩ (φ i (1 -φ j ))Π k (ξ, t) + 1 N ∂ ∂ξ i - ∂ ∂ξ j δ zΩ (ξ i (1 -φ j ) + ξ j φ i )Π k (ξ, t) + 1 2 ∂ ∂ξ i - ∂ ∂ξ j 2 δ zΩ (φ i (1 -φ j ))Π k (ξ, t) + O 1 N 3 2
.

(2.22)

Diffusion in crowded environments

We stress that the variable m i /N can not be replaced by a continuum density because the continuum limit is not appropriate to model the evolution of a single tracer. The following step is to define a new probability distribution Π k (ξ, τ ), which is a function of the vector ξ and of the scaled time τ = t N Ω . In terms of Π k the left hand side of (2.21) becomes

∂P k ∂t = - 1 √ N Ω Ω i=1 ∂Π k ∂ξ i φi + 1 N Ω ∂Π k ∂t .
The leading order contribution in ( 1 √ N ) gives:

-

1 Ω Ω i=1 ∂Π k ∂ξ i φi = δ zΩ Ω i=1 j∈{i-1,i+1} φ i (1 -φ j )( ∂Π k ∂ξ i - ∂Π k ∂ξ j ) (2.23) which yields Ω i=1 ∂Π k ∂ξ i φi = δ z Ω i=1 - ∂Π k ∂ξ i (2φ i -φ i-1 -φ i+1 ), (2.24)
and finally allows one to obtain:

φi = δ 2 ∆φ i . (2.25)
Here ∆ stands for the discrete Laplacian operator defined as:

∆φ i = 2 z j∈i (φ j -φ i ),
with the sum j∈i running over the sites j, which are nearest-neighbors of site i. By taking the size of the patches to zero, one recovers the standard diffusion equation for the density φ of the crowders, in agreement with the result obtained with the other approaches. Taking the following identities:

δ z Ω i=1 j∈{i-1,i+1} ( ∂ ∂ξ i - ∂ ∂ξ j )(ξ i (1 -φ j ) -ξ j φ i )Π k = δ z Ω i=1 ( ∂ ∂ξ i - ∂ ∂ξ i-1 )(ξ i (1 -φ i-1 ) -ξ i-1 φ i )Π k + δ z Ω i=1 ( ∂ ∂ξ i - ∂ ∂ξ i+1 )(ξ i (1 -φ i+1 ) -ξ i+1 φ i )Π k = δ z Ω i=1 ∂ ∂ξ i (ξ i (1 -φ i-1 ) -ξ i-1 φ i )Π k + (ξ i (1 -φ i+1 ) -ξ i+1 φ i )Π k - δ z Ω i=1 ∂ ∂ξ i (ξ i+1 (1 -φ i ) -ξ i φ i+1 )Π k - δ z Ω i=1 ∂ ∂ξ i (ξ i-1 (1 -φ i ) -ξ i φ i-1 )Π k = δ 2 Ω i=1 -∆ξ i Π k and δ z Ω i=1 j∈{i-1,i+1} ∂ ∂ξ i - ∂ ∂ξ j 2 (φ i (1 -φ j ))Π k = δ z Ω i=1 ∂ 2 ∂ξ 2 i φ i (1 -φ i-1 )Π k + ∂ 2 ∂ξ 2 i φ i (1 -φ i+1 )Π k ∂ 2 ∂ξ 2 i+1 φ i (1 -φ i+1 )Π k + ∂ 2 ∂ξ 2 i-1 φ i (1 -φ i-1 )Π k - 2∂ 2 ∂ξ i ∂ξ i-1 φ i (1 -φ i-1 )Π k - 2∂ 2 ∂ξ i ∂ξ i+1 φ i (1 -φ i+1 )Π k = δ z Ω i=1 ∂ 2 δξ 2 i (2φ i + φ i-1 + φ i+1 -2φ i (φ i+1 + φ i-1 ))Π k + ∂ 2 ∂ξ i ∂ξ i-1 -2φ i (1 -φ i-1 ) Π k + ∂ 2 ∂ξ i ∂ξ i+1 -2φ i (1 -φ i+1 ) Π k
at the next-to-leading corrections in eq. (2.21) eventually leads to:

∂Π k ∂t = δ 2 Ω i=1 ∂ ∂ξ i -∆ξ i Π k + δ 2z Ω i=1 i+1 i=i-1 ∂ ∂ξ i ∂ ∂ξ j B i,j Π k + δ z (1 -φ k )Π k-1 -(2 -φ k+1 -φ k-1 )Π k + (1 -φ k )Π k+1 .
Here B is the diffusion matrix, whose entries are

B i,i = 2φ i + φ i-1 + φ i+1 -2φ i (φ i+1 + φ i-1 ) B i,i-1 = (-2φ i (1 -φ i-1 )) B i,i+1 = (-2φ i (1 -φ i+1 )).
(2.26)

To provide a mean-field description we have to consider the probability function of the tagged agent averaged over the fluctuations of the N -particles. This amounts to consider the quantity:

ρ k (t) = Π k dξ
that obeys the following equation:

∂ρ k ∂t = δ z (1 -φ k )ρ k-1 -(2 -φ k+1 -φ k-1 )ρ k + (1 -φ k )ρ k+1 = δ 2 (∆ρ k -φ k ∆ρ k + ρ k ∆φ k ).
The last expression involves the discrete laplacian ∆ defined above. In the continuum limit, if one considers a straightforward generalization to higher dimensions of the previous calculation, one gets :

∂ρ(r, t) ∂t = D ρ (1 -φ(r, t))∇ 2 ρ(r, t) + D ρ ρ(r, t)∇ 2 φ(r, t),
where D ρ is the diffusion coefficient of the tagged particle. The obtained non-linear equation for ρ can be cast as a Fokker-Plank equation:

∂ρ(r, t) ∂t = ∇ 2 D(1 -φ(r, t))ρ(r, t) + 2D∇ ρ(r, t)∇φ(r, t) .
Hence, by neglecting the role of fluctuations, i.e. by operating in the mean-field limit, even in higher dimensions one finds a nonlinear partial differential equation for the density of the tagged species coupled to a standard diffusion equation for the background density:

           ∂φ(r, t) ∂t = D∇ 2 φ(r, t) ∂ρ(r, t) ∂t = ∇ 2 D(1 -φ(r, t))ρ(r, t) + 2D∇ ρ(r, t)∇φ(r, t) .
(2.27)

This system constitutes the generalization of model (2.8) to higher dimensions. It is worth emphasizing that the second equation in (2.27) can be also cast in form of a continuity equation:

∂ρ(r, t) ∂t = ∇ • (D (1 -φ(r, t)) ∇ρ(r, t) + Dρ(r, t)∇φ(r, t)) .
The last procedure that we described is particularly interesting despite its higher complexity for two main reasons. First of all, it allows to extend the description to higher dimensions with respect to the 1D setting where the traditional SEPs are modeled. This is a matter of great importance since to describe the molecular diffusion inside the cell, and in general inside real media, a three dimensional approach is needed. The second reason that contributes to make the Van Kampen's approach interesting is the potential further analysis. More specifically, starting from the microscopic model, in this framework it is possible to go beyond the mean field analysis and monitor the evolution of the fluctuations, which arise due to the finite size effects. However, this analysis lies outside the purpose of this work, that aims at characterizing the effect of crowding on the macroscopic quantities. To validate the soundness of the mean-field approach in two dimension for the aforementioned setting, we will provide hereafter a brief comparison between the mean-field prediction in 2D with the stochastic model.

On the validity of the mean-field approximation: comparing stochastic and mean-field simulations in 2D

To test the adequacy of the proposed mean-field model we have performed a campaign of numerical simulations, with reference to the 2D setting. More specifically we have implemented a Monte Carlo scheme to solve the stochastic process under scrutiny and so trace the evolution of the tagged particle in time. At each time iteration, all crowders and the tagged particle can update their position, moving at random, and with equal probability, in one of the four allowed directions, provided the selected target site is unoccupied. If the destination site is occupied, the move is rejected and the particles keep their original positions. The order of selection of the particles is, at each iteration, randomized. By averaging over many independent realizations, one can reconstruct the normalized histogram of the position visited by the tagged particles at a given time t and compare it with the density profile ρ obtained upon integration of the mean-field system (2.27). To carry out the numerical integration of the above partial differential equations we assumed a forward difference approximation in time and replaced the spatial derivatives by centered approximations. The result of the comparison is reported in figures 2.4(a) and 2.4(b), for two choices of the initial conditions. In figure 2.4(a), the tagged particle is initially positioned in the middle of a two dimensional waterbag, filled with crowders, with average density equal to φ 0 . In figure 2.4(a), the tagged particle is instead positioned, at time t = 0, in the center of an empty region, a square of assigned size. The crowders are instead assumed to occupy an adjacent domain with average uniform density φ 0 . In both cases, the agreement between stochastic and mean field simulations is satisfying.

Evolution of the mean square displacement: sub-and superdiffusive transients

Once the equations (2.8) describing the evolution of the density of the tagged species and the bulk population have been recovered, we shall use the results to identify the effect of the crowding on the global kinetics. More specifically, to understand how it impacts the motion of the tagged species, we study the time evolution of the mean square displacement (MSD) µ 2 (t) of the recognizable particle, defined as:

µ 2 (t) = ρ(x, t)|x -x | 2 d n x (2.28)
where x = ρ(x, t)x d n x, and the superscript n denotes the space dimension.

It is well known that the MSD scales linearly with time for unobstructed diffusion, while a sub-linear growth of the MSD is often interpreted as a direct manifestation of the microscopic competition for available space in crowded media. As we shall prove in the following, this is an overly simplistic picture, as more complex scenarios can easily be obtained by direct integration of Eqs. (2.8), where nonlinear MSDs emerge only as transient regimes. We are particularly interested in a specific class of initial condition, symmetric in the domain of definition, so that x = 0.

Let us first illustrate the one-dimensional case. At time t = 0, the tagged species is localized at the origin, while the crowders populate a compact domain also centered at the origin. In formulae, ρ(x, 0) = δ(x), where δ(•) is Dirac delta, and φ(x, 0) = φ 0 [θ(x + x 0 )θ(xx 0 )], with φ 0 ∈ [0, 1] gauging the crowding strength (see inset in the upper panel of Fig. 2.5). From here on, as a further simplification, we assume D ρ = D φ = D. Fig. 2.5 shows the rescaled MSD of the tagged species as a function of time as obtained by integrating Eqs. (2.8) numerically. At short times, the tagged species is immersed in the almost uniform sea of surrounding crowders. Since φ is approximately constant, the tagged particles diffuse normally with an effective diffusion coefficient equal to D(1φ 0 ). In the long-time limit, the crowders are evenly spread over the one-dimensional support (which we imagine closed but very large so as to neglect boundary effects). Consequently, the density φ is small and its contribution can be neglected in the Smoluchowski equation for the evolution of ρ (zero Diffusion in crowded environments It is remarkable, and to some extent counter-intuitive, that a super-diffusive transient is found in a model accounting for crowding in the absence of driving. In fact, excluded-volume interactions among diffusing agents are customarily believed to lead to slower-than-diffusive spread of concentrations. We observe that the time duration of the super-diffusive transient increases quadratically with x 0 , the width of the initial density of crowders. Of course, the origin of the observed dynamics can be traced back to the force term in the current of tagged particles. The effective force F (x, t) = -∂φ(x, t)/∂x induced by the crowders acts as a systematic bias in the evolution of the density ρ. Initially, ρ evolves freely, as φ(x, t) ≃ φ 0 for all values of x where ρ is non-zero. Then, after a time of the order of τ ∝ x 2 0 , the support of ρ extends to a domain where it is no longer possible to assume φ(x, t) constant. In particular, ∂φ(x, t)/∂x < 0 for x > 0 and ∂φ(x, t)/∂x > 0 when x < 0, which implies F (x, t) > 0, hence a forceinduced boost over the osmotic current. The mean-field force, which stems from the microscopic competition for space between crowders and tagged particles, pulls the distribution ρ away from the origin, stretching the right (left) tail towards the direction of positive (negative) x. This leads to the super-diffusive transient shown in the upper figure of Fig. 2.5.

A dual situation can be imagined yielding a sub-diffusive transient. To this end, let us consider the crowders to be initially distributed uniformly in a (large) one-dimensional domain of size 2L. At time t = 0, the crowders that populate a segment of width 2x 0 , centered around the origin, are removed from the system. This amounts to considering the initial distribution φ(x, 0) = φ 0 [1 -θ(x + x 0 ) + θ(xx 0 )] (see inset in the bottom panel of Fig. 2.5). At short times, the diffusion of tagged particles is not affected by the crowders. The rescaled MSD µ 2 /2Dt is hence approximately equal to one and stays constant over a finite time window of order x 2 0 /D. At long times, the crowders will have approximately relaxed to the uniform concentration φ L = φ 0 (1 -x 0 /L). Hence, the tagged particles will find themselves diffusing in a uniform medium with a reduced diffusion coefficient D(1 -φ L ). The bottom panel of Fig. 2.5 confirms our reasoning, as the rescaled MSD µ 2 /2Dt is seen to decrease monotonously, interpolating between the initial plateau µ 2 /2Dt = 1 and the final value µ 2 /2Dt = (1 -φ L ) < 1. In this case, one thus observes a sub-diffusive crossover. In fact, in this case ∂φ(x, t)/∂x > 0, for x > 0, which implies F (x, t) < 0, i.e. an effective force that opposes the osmotic thrust to delocalization.

Let us now turn to considering the spreading of tagged particles in two dimension. To this end, we consider a straightforward generalization of the initial conditions discussed above, as exemplified by the cartoons reported in Figure 2.6. Sub and super-diffusive transients are again observed depending on the initial conditions, in stringent analogy with what observed in one dimension. Along the same lines, the observed behavior can be rationalized in terms of the effective force on the tagged species caused by the crowders. In this case, the total radial current is J r = -D{∇ r [(1 -φ)ρ] + 2ρ∇ r φ}. Again, we see that an effective force in the radial direction F = 2∇ r φ arises when the tagged particles diffuse into regions of changing density of crowders, yielding a current boost (super-diffusive crossover) or a drop (sub-diffusive crossover) depending on the initial conditions.

The situations described above can be easily recreated in laboratory experiments, by initially confining the particles, including those whose evolution is to be tracked, within a finite portion of the available space.

In summary, our results prove that both super-diffusion and sub-diffusion transients can occur as a result of crowding in one and higher dimensions in the absence of driving, depending on the initial conditions. These findings strongly warn against the simplistic identification of crowding with anomalous transport tout court.

Diffusion in crowded environments

The persistent random walk with exclusion and its continuum limit

In this section we want to introduce a different rule in the microscopic description of a random walk, prescribing that the agents, in addition to being subject to the excluded volume effects, also have a bias to keep hopping in the same direction as they did in the past. A random walk in which the walker has probability α of continuing in the same direction as the previous step is known as persistent random walk (PRW). While the continuum limit of standard RWs is the diffusion equation, which yields an infinite propagation velocity, the continuum limit of the PRW is the so-called telegraph equation, which displays a transition from ballistic to diffusive transport at a characteristic time. The PRW and its connection with the telegraph equation were first studied by Goldstein in 1951 [START_REF] Goldstein | On diffusion by discontinuous movements, and on the telegraph equation[END_REF]. An interesting discussion of many applications of the telegraph equation can be found in a recent review by G. Weiss [START_REF] George | Some applications of persistent random walks and the telegrapher's equation[END_REF].

Working again in the same framework of SEPs models, we generalize the concept of persistent random walk to the case of interest where exclusion effects are to be accounted for. In the classical microscopic formulation of the telegraph equation, individual walkers are assumed to jump towards neighboring sites, with constant probability. One could imagine to modify the PRW by introducing an explicit constraint in the probability of jumps that weights the occupancy of the target sites. In doing so, we will obtain a generalized telegraph equation which includes non linear terms, reflecting the microscopic competition for the available spatial resources. It will appear how, while the nonlinear terms prove negligible in the diluted limit, when working at high densities excluded-volume corrections do matter. We will substantiate this claim both analytically and numerically.

From here on we will refer to the generalized model here introduced as to a persistent simple exclusion process (PSEP).

Mathematical model for the persistent random walk with exclusion

Let us consider a number N walkers on a one-dimensional lattice with spacing d and length L. According to the definition of persistent random walk [START_REF] Goldstein | On diffusion by discontinuous movements, and on the telegraph equation[END_REF], at regular intervals ∆t a walker can jump in the same direction as it did at the previous step with probability p or invert its direction with probability q. We take q = 1-p, which amounts to assuming that there is no leakage [START_REF] Goldstein | On diffusion by discontinuous movements, and on the telegraph equation[END_REF] in the system. Let us denote with a i (n) the probability that a walker is at site i at time n∆t having been at site i -1 at time (n -1)∆t (right-bound flow) and with b i (n) the probability that a walker is at site i at time n∆t having been at site i + 1 at time (n -1)∆t (left-bound flow).

If walkers are invisible to each other, the following relations hold

a i (n) = p a i-1 (n -1) + q b i-1 (n -1) (2.29) b i (n) = p b i+1 (n -1) + q a i+1 (n -1).
(2.30)

The above equations describe a discrete stochastic process. The continuum limit can be obtained by introducing the continuous probability density

P (x, t) = P i (n) ≡ a i (n) + b i (n) ,
where . . . denotes average over the trajectory of many agents, obtained by letting d → 0, ∆t → 0, q → 0. Working in this framework, it is known that one gets the telegraph equation [82]

∂ 2 P ∂t 2 + 2r ∂P ∂t = v 2 ∂ 2 P ∂ 2 x (2.31) with lim d,∆t→0 d ∆t = v lim q,∆t→0 q ∆t = r. (2.32)
In this section we wish to study how a persistent random walk is modified by enforcing the constraint that prevents two walkers to occupy the same site at the same time. That is, if the probability to jump to a given site is gauged by the current occupancy of that site. Along the same line of reasoning of SEPs and ASEPs, we modify Eqs. (2.29) and (2.30) in the following way:

a i (n) -a i (n -1) = [p a i-1 (n -1) + q b i-1 (n -1)][1 -P i (n)] (2.33) -a i (n -1){p[1 -P i+1 (n -1)] + q[1 -P i-1 (n -1)]} b i (n) -b i (n -1) = [p b i+1 (n -1) + q a i+1 (n -1)][1 -P i (n)] (2.34) -b i (n -1){p[1 -P i-1 (n -1)] + q[1 -P i+1 (n -1)]} where P i (n) = a i (n) + b i (n).
Again, the idea is to gauge jump probabilities by the occupancy of the target sites. For example, the first term in the right hand side of eq.(2.33) states that a net increase of the probability at site i associated with the right-bound flow is only possible with a transition rate proportional to the amount of free room at site i, i.e. (1 -P i ). If P i = 1, no further increase of a i nor of b i is possible.

In order to take the continuum limit, we first divide Eqs. (2.33) and (2.35) by ∆t and substitute q = 1p. Then, recalling the definitions (2.32), we get

∂a ∂t + v ∂ ∂x [a(1 -P )] = -rJ(1 -P ) ∂b ∂t -v ∂ ∂x [b(1 -P )] = rJ(1 -P ) (2.35)
where P (x, t) = a(x, t) + b(x, t) and J(x, t) = a(x, t)b(x, t). For the sake of the argument, let us consider the propagation of pulses in a fluid, i.e. traveling density fluctuations. Eqs. (2.35) contain the single-particle probability field P , which is a number between zero and one. The value P = 1 should correspond to the maximum density allowed in the system. Thus, more physical equations can be obtained by introducing the agent densities ρ(x, t) ≡ ρ M P (x, t) (2.36)

J (x, t) ≡ ρ M J(x, t) (2.37) ρ + (x, t) = ρ M a(x, t) (2.38) ρ -(x, t) = ρ M b(x, t) (2.39)
where ρ M is the maximum allowed density, which in principle could be regarded as a parameter of the model. If we imagine that the agents have a finite size σ, i.e. we Diffusion in crowded environments regard them as hard rods, one simply has ρ M = 1/σ. Introducing the density ρ M , eqs.2.35 become

∂ρ + ∂t + v ∂ ∂x ρ + 1 - ρ ρ M = -rJ 1 - ρ ρ M ∂ρ - ∂t -v ∂ ∂x ρ -1 - ρ ρ M = rJ 1 - ρ ρ M .
(2.40)

A system of equations for the densities ρ(x, t) and J (x, t) can be obtained by adding and subtracting the two equations 2.35

∂ρ ∂t + v ∂ ∂x J 1 - ρ ρ M = 0 ∂J ∂t + v ∂ ∂x ρ 1 - ρ ρ M = -2rJ 1 - ρ ρ M (2.41)
As a general remark, we see that the microscopic exclusion constraint results in the appearance of nonlinear terms. The standard evolution of the PRW leading to the telegraph equation is obtained in the dilute limit ρ ≪ ρ M . Conversely, we may consider the full system (2.41) when describing transport in a crowded medium. The nonlinear equations here embody the microscopic excluded-volume constraint that impact strongly the evolution at high densities.

2.3.2

Mean square displacement in the limits t → 0 and t → ∞

We turn now to analyzing how the excluded-volume constraint affects the propagation of an initially localized pulse. It is well known that the PRW displays a transition from ballistic to diffusive transport, as exemplified by the mean square displacement (MSD),

µ 2 (t) ≡ 1 N x 2 (t) ρ -x(t) 2 ρ (2.42) with x m (t) ρ = x m ρ(x, t)dx and N = ρ(x, t)dx.
As it is customarily done, we shall here restrict to a class of symmetric initial pulses, namely such that ρ(x, t = 0) = ρ(-x, t = 0) and J (x, t = 0) = 0. This amounts to considering the initial distribution of the right-headed agents equal to that of the left-headed ones. In this case, it is straightforward to show that x(t) = 0 ∀ t.

For the PRW one has

µ 2 (t) -µ 2 (0) = v 2 2r 2 (2rt -1 + e -2rt ) ≃    v 2 t 2 for t ≪ 1/2r v 2 r t for t ≫ 1/2r
When excluded-volume effects are important, it appears impossible to obtain a closed expression for µ 2 (t). However, one can still capture the asymptotics in a particular symmetric regime and in the absence of boundary terms. Let us consider the Taylor expansion of µ 2 (t)

µ 2 (t) = µ 2 (0) + µ ′ 2 (0) t + 1 2 µ ′′ 2 (0) t 2 + O(t 3 ) (2.

43)

In order to evaluate the coefficients of the expansion, let us multiply the first equation of eq. (2.41) by x 2 and the second one by x and integrate. Integrating by parts and assuming that boundary terms vanish, we obtain

d dt x 2 ρ -2v x J + 2v ρ M x J ρ = 0 d dt x J -v ρ 1 - ρ ρ M dx + 2r x J - 2r ρ M x J ρ = 0 (2.44)
where . . . J and . . . J ρ denote averages with respect to the corresponding (products of) densities.

Since the initial condition is symmetric, the first equation of (2.44) shows that µ ′ 2 (0) = 0. Differentiating the same equation with respect to time, recalling equations 2.41 and integrating by parts eventually leads to

d 2 dt 2 x 2 ρ = 2v d dt x J - 2v ρ M d dt x J ρ = 2v v ρ 1 - ρ ρ M dx -2r x J + 2r ρ M x J ρ - 2v ρ M v ρ 2 1 - ρ ρ M dx + v xρ 1 - ρ ρ M ∂ρ ∂x dx (2.45)
Evaluating the previous expression at t = 0 makes the terms involving J disappear. Thus

d 2 dt 2 x 2 ρ t=0 = 2v 2 ρ 1 - ρ ρ M 2 dx t=0 - 2v 2 ρ M xρ 1 - ρ ρ M ∂ρ ∂x dx t=0 (2.46)
which yields:

µ 2 (t) ≈ µ 2 (0) + v 2 e t 2 (2.47)
with

v e = v √ N ρ 1 - ρ ρ M 2 dx t=0 - 1 ρ M xρ 1 - ρ ρ M ∂ρ ∂x dx t=0 1/2 (2.48)
where we have used the fact that the quantity N = ρ(x, t) dx is constant if the boundary terms can be neglected. This is the case corresponding to a broad initial conditions, such as in the propagation of pulses initially displaying a finite support. Fig. 2.7 shows the time evolution of µ 2 obtained by integrating numerically the system (2.44) with a forward-difference approximation in time and replacing the spatial derivatives by centered Euler approximations. The initial conditions are generalized Gaussian pulses of the type ρ(x, t = 0) = φ e -x 2β /2σ 2 .

(2.49)

This choice allows us to investigate the propagation of a pulse whose shape varies continuously from Gaussian (β = 1) to a sharp step function β ≫ 1, while the parameter φ < 1 gauges the crowding. It is not difficult to compute v e analytically The solid line refers to the PRW and vanishes as t -1 . The dashed line is an inversepower law with exponent 0.4, to be used as a guide for the eye.

from eq. (2.48) as a function of φ for a pulse of the kind (2.49). After straightforward calculations, one gets

v e (φ) = v 1 - 3 φ 2 γ + 2 φ 2 3 γ 1/2 (2.50)
where γ = 1 + 1/2β. It is clear from the figure that the approximation (2.50) captures to an excellent extent the initial ballistic stage. Furthermore, the numerical integration of eqs. 2.44 shows that asymptotically the propagation becomes diffusive, with the same diffusion coefficient v 2 /2r as the PRW. This is to be expected as ρ → 0 as t → ∞ and therefore the excluded-volume constraints (that is, the nonlinear terms) become negligible. Nevertheless, the inset in fig. 2.7 clearly shows that the approach to the diffusive regime is considerably slowed down as a result of crowding, and the more the greater the excluded-volume constraint.

Our analysis shows that the initial behavior of the mean square displacement is qualitatively the same as in the PRW, i.e. the propagation is ballistic. The effect of crowding is to decrease the velocity that characterizes the initial stage of the evolution. In Fig. 2.8 the effective velocity v e , normalized to the diluted limit v, is plotted as function of the level of crowding and for different choices of the parameter β. Interestingly, v e depends on the shape of the initially localized density pulse. The ballistic spreading of a super-Gaussian pulse, nearly a sharp step, proceeds with a considerably lower speed as compared to the spreading of a pure Gaussian pulse (see again eq. (2.48)). At large times, the mean square displacement for the telegraph equation becomes linear. Our goal is to understand whether the excluded-volume effects alter this behavior. Let us suppose that (ρ, J ) tends to a stationary solution for the system (2.41) as t goes to ∞. Then we have:

lim t→∞ ρ(x, t) = ρ ∞ lim t→∞ J (x, t) = 0
where ρ ∞ depends on the size of the domain since ρ(x, t)dx = C, with the constant C fixed by the initial conditions. Indeed

C = lim t→∞ L 2 -L 2 ρ(x, t)dx = ρ ∞ L. (2.51)
The effect of the term due to the crowding ρ 2 dx within the system (2.44) is negligible with respect to ρ(x, t)dx as the domain becomes larger

lim t→∞ L 2 -L 2 ρ 2 (x, t)dx L 2 -L 2 ρ(x, t)dx = C L ≪ 1 if L ≫ 1. (2.52)
Diffusion in crowded environments

Hence the system (2.44) for t ≫ 1 is approximated by:

d dt x 2 ρ -2v[z] = 0 d dt x J -vC + 2r[z] = 0
where z = x J -2 ρ M x J ρ . If we make the assumption that x J tends to some constant value, we get from the previous system

d dt x 2 ρ = v 2 C r (2.53)
resulting in the same expression for the mean square displacement in the ∞-limit that we find in the diluted case

µ 2 (t) ∼ v 2 r t + µ 2 (0) if t ≫ 1/2r. (2.54)

Discussion on the emergence of transient regimes

In last two sections we have analyzed the excluded-volume effects on two examples of random walks governed by different microscopic rules. Working with a tagged particle, we have shown that the excluded-volume effect can emerge in the macroscopic picture as a transient regime, either sub-diffusive or super-diffusive, depending on the specific initial conditions. When the crowders are uniformly dispersed around an isolated patch where the tagged species is initially confined, a sub-diffusive scaling for the mean square displacement of the density of the tagged species is observed. This crossover regime persists within a finite, possibly very long time window. On the contrary, if the tagged agents are trapped inside a uniform patch of crowders confined in a much larger, otherwise empty container, the excluded-volume interactions produce an effective force term in the current of tagged particles, that accelerates their spread with respect to the osmotic current, yielding a super-diffusive transient. The effect of the crowding is thus temporary, as standard diffusion is recovered in the long time limit. A similar temporary effect is observed in the analysis of the persistent random walk with exclusion. The continuous limit of the persistent random walk with exclusions undergoes a transition from ballistic to diffusive regimes, showing the same transition as the corresponding persistent random walk which holds under diluted conditions. Thus, from a qualitative point of view, both the SEP and the PSEP models yield the same behavior. Nevertheless, by evaluating the parameters, we pointed out that the crowding causes the velocity in the ballistic phase to decrease. The effective velocity depends on the shape of the initially localized density pulse: the ballistic spreading of a sharp step function proceeds with a considerably lower velocity, as compared to that associated with a pure Gaussian pulse. Asymptotically, the propagation becomes diffusive, and it displays the same diffusion coefficient as the PRW, making the excluded-volume effect disappear in the long limit. We stress that crowding affects the macroscopic picture by inducing transitions between nonanomalous (standard) evolution regimes. These transitions can be mistaken for more general hallmarks if one restricts himself to insufficiently long time series.

Excluded-volume effects in non-homogeneous media

In the previous description we have described stochastic processes with homogeneous jump probabilities. More specifically, by taking q j = q in all the described processes, we assumed that diffusion did not depend on the agent's position within the domain. This amounts to model diffusion of particles in homogeneous environments, and eventually leads to a constant diffusion coefficient in the macroscopic picture. We will now move on to analyze diffusion processes inside inhomogeneus media. This, at the microscopic level, amounts to considering stochastic processes where the jumps are controlled by site-dependent hopping rates. We will show how, enforcing these rules, one can obtain three different nonlinear advection-diffusion equations in the continuum limit, depending on the particular rule that is used to embody site inequivalence at the microscopic level. Moreover we will see that the macroscopic description changes whether we take into account or not the excluded volume effect.

A stochastic jump process on a 1D lattice with inequivalent sites in the presence of a field and with excluded volume effect describes the most general SEP in one dimension:

n i (k + 1) -n i (k) = z + i-1 n i-1 (k)[1 -n i (k)] + z - i+1 n i+1 (k)[1 -n i (k)] -z + i n i (k)[1 -n i+1 (k)] -z - i n i (k)[1 -n i-1 (k)].
(2.55) Similarly to the previous processes, eq. (2.55) is to be regarded as the update rule for a Monte Carlo process, where n i (k) is the occupancy of site i at time t = k∆t, which can be either zero or one. The variables z ± i take the value 0 or 1 according to a random number ξ i which has a uniform distribution between 0 and 1. By defining the jump probabilities q ± j (j = i, i ± 1) one can formally write:

z + i-1 = θ(ξ i ) -θ(ξ i -q + i-1 ) z - i+1 = θ(ξ i -q + i-1 ) -θ(ξ i -q + i-1 -q - i+1 ) z + i = θ(ξ i -q + i-1 -q - i+1 ) -θ(ξ i -q + i-1 -q - i+1 -q + i ) z - i = θ(ξ i -q + i-1 -q - i+1 -q + i ) -θ(ξ i -1) (2.56)
where θ(•) stands for the Heaviside step function and where we are assuming that

q + i-1 + q - i+1 + q + i + q - i = 1.
The ordering of appearance of the q ± j in the above expressions is arbitrary. Equations (2.56) entail that z ± j = q ± j , where • indicates an average over many values of ξ i , for a given configuration {n i }. The above process is fully determined by the fields q ± i , specifying the probability of jumping from site i to site i + 1 (q + i ) or from site i to site i -1 (q - i ) in a time interval ∆t. A (discrete-time) master equation for the above SEP can be obtained by averaging over many Monte Carlo cycles performed according to rule (2.55)

P i (k + 1) -P i (k) = q + i-1 [P i-1 (k) -P i,i-1 (k)] + q - i+1 [P i+1 (k) -P i,i+1 (k)] -q + i [P i (k) -P i,i+1 (k)] -q - i [P i (k) -P i,i-1 (k)] (2.
57) where we have defined the one-body and two-body site occupancy probabilities

P i (k) = n i (k) (2.58) P i,i±1 (k) = n i (k)n i±1 (k) (2.

59)

Here • denotes averages performed over many independent Monte Carlo cycles performed until time k∆t starting from the same initial condition. We emphasize that the same equation has been derived through a slightly different procedure by Richards in 1977 [START_REF] Peter | Theory of one-dimensional hopping conductivity and diffusion[END_REF].

A pondering pause is required at this point before carrying out the continuum limit. In fact, a moment's thought is enough to realize that there are (at least) three ways one can enforce quenched disorder, corresponding to spatially varying hopping rates, in a jump process with the aim of modeling propagation in a inhomogenous medium. The master equation (2.57) reflects only one of the possible choices.

Let us consider the jump from site i to site i + 1. The probability of an agent taking an i → i + 1 leap can be equally well taken as (1) proportional to q + i , (2) proportional to q + i+1 or (3) proportional to (q + i + q + i+1 )/2. Let a be the lattice spacing and let us define a reversal probability ǫ i , such that

q + i = Q i q - i = Q i -ǫ i (2.60)
The condition (2.60) (with ǫ i > 0) amounts to considering a field introducing a bias in the positive x direction. In order to take the continuum limit lim a,∆t→0 P i (k) = P (x, t), we must require lim a,∆t→0

Q i a 2 ∆t = D(x) (2.61) lim a,∆t→0 ǫ i a ∆t = v(x) (2.62)
Eq. (2.61) defines the position-dependent diffusion coefficient, while eq. (2.62) defines the field-induced drift velocity. Note that we are assuming that the reversal probability vanishes linearly with a.

All three cases correspond to the same space-dependent function D(x) in the continuum limit, as prescribed by eq. (2.61). However, as we shall see in the following, depending on whether rule (1), ( 2) or ( 3) is chosen, one is led to totally different advection-diffusion equations in the continuum limit.

Processes without exclusion

In order to illustrate this subtle point, let us start with jump processes in the presence of quenched disorder but with no exclusion constraints on the allowed moves. Following the same reasoning on the correlation probability, it is not difficult to realize that the three possible choices (1), ( 2) and (3) referred above lead to the following master equations

P i (k + 1) -P i (k) = q - i+1 P i+1 (k) + q + i-1 P i-1 (k) -(q + i + q - i )P i (k) (2.63) P i (k + 1) -P i (k) = q - i P i+1 (k) + q + i P i-1 (k) -(q + i+1 + q - i-1 )P i (k) (2.64) P i (k + 1) -P i (k) = q - i+1 + q - i 2 P i+1 (k) + q + i-1 + q + i 2 P i-1 (k) -P i (k) q + i+1 + q + i 2 + q - i-1 + q - i 2 (2.65)
In the first case, eq. (2.63) (already obtained above), the rate for a given jump depends on the starting site, while the second, eq. (2.64), is the opposite, the rate depends on the index of the target site. The third case, eq. (2.65) is an intermediate, symmetric situation where jump rates are associated with links rather than with nodes. In the continuum limit a, ∆t → 0 one has

Q i±1 ≈ Q(x) ± ∂Q ∂x a + 1 2 ∂ 2 Q ∂x 2 a 2 + . . . ǫ i±1 ≈ ǫ(x) ± ∂ǫ ∂x a + . . . P i±1 (k) ≈ P (x, t) ± ∂P ∂x a + 1 2
∂ 2 P ∂x 2 a 2 + . . . ( 2
∂P (x, t) ∂t = - ∂J 1 (x, t) ∂x (2.67) ∂P (x, t) ∂t = - ∂J 2 (x, t) ∂x (2.68) ∂P (x, t) ∂t = - ∂J(x, t) ∂x (2.69)
where

J 1 (x, t) = -D(x) ∂P (x, t) ∂x + v(x) - ∂D(x) ∂x P (x, t) (2.70) J 2 (x, t) = -D(x) ∂P (x, t) ∂x + v(x) + ∂D(x) ∂x P (x, t) (2.71) J(x, t) = -D(x) ∂P (x, t) ∂x + v(x)P (x, t). (2.72)
We see that the stochastic processes (2.63), (2.64) and (2.65) correspond to diffusion with drift (or, equivalently, in an external potential). In the cases of the two processes (2.63) and (2.64) the drift velocity comprises two contributions: the difference between the right-bound and left-bound jump rate fields and a contribution arising from the spatial variation of the diffusion coefficient. Interestingly, the latter term has the opposite sign depending on whether jumps at the microscopic level are controlled by the rates evaluated at the start or at the target sites. Conversely, considering symmetrized jump rates does not result in the appearance of such additional term in the drift velocity.

Enforcing the excluded-volume constraint

We can now come back to our original aim, i.e. taking the continuum limit of the master equation (2.55). It is now clear that, if we want to consider an inhomogeneous medium, we must not restrict to the prescription leading to eq. (2.55), but we must also consider the other two kinds of processes described above in the absence of exclusion.

Extending the reasoning that led us to eq. (2.55) and assuming the customarily mean-field (MF) factorization from the microscopic stochastic process described by

P i,i±1 (k) ≡ n i (k)n i±1 (k) = n i (k) n i±1 (k) = P i (k)P i±1 (k), (2.73)
we obtain:

P i (k + 1) -P i (k) = [q - i+1 P i+1 (k) + q + i-1 P i-1 (k)][1 -P i (k)] -P i (k){q + i [1 -P i+1 (k)] + q - i [1 -P i-1 (k)]} (2.74) P i (k + 1) -P i (k) = [q - i P i+1 (k) + q + i P i-1 (k)][1 -P i (k)] -P i (k){q + i+1 [1 -P i+1 (k)] + q - i-1 [1 -P i-1 (k)]} (2.75) P i (k + 1) -P i (k) = q - i+1 + q - i 2 P i+1 (k)[1 -P i (k)] + q + i-1 + q + i 2 P i-1 (k)[1 -P i (k)] - q + i+1 + q + i 2 P i (k)[1 -P i+1 (k)] - q - i-1 + q - i 2 P i (k)[1 -P i-1 (k)]. (2.76)
The continuum limit a, ∆t → 0 of the above master equations is readily obtained by introducing as above the Taylor expansions of q(x) and P (x, t). In this way we find:

∂P (x, t) ∂t = - ∂J 1 (x, t) ∂x (2.77) ∂P (x, t) ∂t = - ∂J 2 (x, t) ∂x (2.78) ∂P (x, t) ∂t = - ∂J 3 (x, t) ∂x (2.79)
where

J 1 (x, t) = -[1 -P (x, t)] ∂ ∂x [D(x)P (x, t)] -D(x)P (x, t) ∂P (x, t) ∂x + v(x)P (x, t)[1 -P (x, t)]
(2.80)

J 2 (x, t) = P (x, t) ∂ ∂x {D(x)[1 -P (x, t)]} -D(x)[1 -P (x, t)] ∂P (x, t) ∂x + v(x)P (x, t)[1 -P (x, t)]
(2.81)

J 3 (x, t) = -D(x) ∂P (x, t) ∂x + v(x)P (x, t)[1 -P (x, t)]. (2.82)
We see that even in the case of excluded-volume interactions, the mean-field equations can be cast in the form of continuity equations with suitably defined currents given by eqs. (2.80), (2.81) and (2.82). Moreover, as for the symmetric exclusion process without quenched disorder, the microscopic exclusion constraint disappears in taking the continuum limit of the master equation (2.76) for q + i = q - i , which yields a zero-advection-diffusion equation, eq. (2.79), identical to its counterpart with no exclusion, eq. (2.69).

Eqs. (2.77), (2.78) and (2.79) are nonlinear advection-diffusion equation, appropriate for describing the continuum limit of a microscopic exclusion process occurring on a lattice of inequivalent sites in the presence of a field. It is interesting to note that in the case of equivalent sites, which translates to a constant diffusion coefficient, the diffusive parts become linear, i.e. the microscopic exclusion rule is lost in the diffusive part. In the case of zero field, one then simply recovers the ordinary diffusion equation for the three jump processes, which, as it is widely known, can be derived from a microscopic jump process with no exclusion rules. This curious observation has been first reported by Huber [START_REF] Huber | Particle kinetics on one-dimensional lattices with inequivalent sites[END_REF]. If both the diffusion coefficient and the drift velocity are constant, eqs. (2.77), (2.78) and (2.79) all reduce to

∂P ∂t = D∇ 2 P -v ∂ ∂x [P (1 -P )] (2.83)
an equation already obtained recently in Ref. [START_REF] Simpson | Pathlines in exclusion processes[END_REF].

Eqs. (2.77), (2.78) and (2.79) contain the single-particle probability field P (x, t), which is a number between zero and one. The value P = 1 should correspond to the maximum density ρ M allowed in the system. Thus, more physical equations can be obtained as in section 2.3 by introducing the agent density

ρ(x, t) ≡ ρ M P (x, t) = φ M v d (σ/2) P (x, t) (2.84)
where

v d (r) = (π 1/2 r) d Γ(1 + d/2) (2.85)
is the volume of a d-dimensional sphere1 of radius r and φ M is the maximum packing fraction for systems of d-dimensional hard spheres,

φ M = 1 (d = 1), φ M = π/ √ 12 ≈ 0.907 (d = 2) and φ M = π/ √ 18 ≈ 0.740 (d = 3) [85]
. With these definitions, and using a more general vector notation, eqs. (2.77), (2.78) and (2.79) become

∂ρ(x, t) ∂t = -∇ • J 1 (x, t) ∂ρ(x, t) ∂t = -∇ • J 2 (x, t) ∂ρ(x, t) ∂t = -∇ • J 3 (x, t) (2.86) with J 1 (x, t) = -1 - ρ ρ M ∇[D(x)ρ(x, t)] -D(x) ρ ρ M ∇ρ(x, t) + v(x)ρ(x, t) 1 - ρ ρ M J 2 (x, t) = ρ(x, t)∇ D(x) 1 - ρ ρ M -D(x) 1 - ρ ρ M ∇ρ(x, t) + v(x)ρ(x, t) 1 - ρ ρ M J 3 (x, t) = -D(x)∇ρ(x, t) + v(x)ρ(x, t) 1 - ρ ρ M .
For a homogeneous medium in the absence of a field these reduce to a simple diffusion equation. Only for inequivalent sites and/or in the presence of a field the microscopic exclusion constraint does survive in the mean-field limit.

As a final remark, we note that analogous models have been derived following similar lines in Ref. [START_REF] Lakatos | Hydrodynamic mean-field solutions of 1D exclusion processes with spatially varying hopping rates[END_REF][START_REF] Kevin | Modelling the movement of interacting cell populations[END_REF]. Nonlinear mean-field equations for exclusion process of this type have been used since the 70s to investigate one-dimensional transport in solids [START_REF] Huber | Particle kinetics on one-dimensional lattices with inequivalent sites[END_REF]. In fact, despite mean-field descriptions for the inhomogeneous ASEP are known to provide imperfect descriptions of certain non-equilibrium observables in one dimension, e.g. the current-density relation and critical exponents [START_REF] Szavits | Disordered exclusion process revisited: some exact results in the low-current regime[END_REF], continuum descriptions can be employed reliably to track the time-evolution of largewavelength density fluctuations [START_REF] Lakatos | Hydrodynamic mean-field solutions of 1D exclusion processes with spatially varying hopping rates[END_REF][START_REF] Galanti | Diffusion of tagged particles in a crowded medium[END_REF][START_REF] Fernando | Nonlinear diffusion and exclusion processes with contact interactions[END_REF][START_REF] Landman | Myopic random walkers and exclusion processes: Single and multispecies[END_REF][START_REF] Penington | Building macroscale models from microscale probabilistic models: A general probabilistic approach for nonlinear diffusion and multispecies phenomena[END_REF][START_REF] Penington | Collective motion of dimers[END_REF][START_REF] Penington | Interacting motile agents: Taking a mean-field approach beyond monomers and nearestneighbor steps[END_REF][START_REF] Simpson | Pathlines in exclusion processes[END_REF].

Extended crowding

In the previous section we described how to derive macroscopic mean-field equations starting from microscopic master equations which account for exclusion effects. In our description we did not take into account the size and the shape of the agents by prescribing that each molecule can occupy just one site and can not jump into a site which is already occupied by another agent. To recover a macroscopic continuous description, we performed the limit for vanishing lattice spacing. This strategy amounts to considering agents of vanishing size in the continuum limit. We term this peculiar situation in the macroscopic world point-like crowding, as opposed to the concept of extended crowding which takes into account the finite dimension of the particles. As we have shown in the previous sections, considerable microscopic information is lost in the continuum limit with point-like agents. The point-like characterization has to be the reason why the mean-field approximation loses the memory of the microscopic exclusion constraint and the diffusion equation is recovered for equivalent sites in the absence of a field. In this spirit we want to investigate whether working in the framework of extended crowding allows one to keep track of the excluded volume effects even in the absence of external forces or inhomogeneities in the domain. Moreover, in many biological contexts, the different shape and dimension of individuals in various populations which are diffusing in the same environment do have an influence on the extent of the collective motion. For these reasons we are interested in microscopic exclusion processes involving agents characterized by a finite size, as opposed to standard SEPs.

Extended crowding models and possible interpretations

The fact that cells are often elongated and rod-shaped inspired the development of exclusion process models using elongated agents in the form of hard rods as the individual units [START_REF] Simpson | Models of collective cell spreading with variable cell aspect ratio: A motivation for degenerate diffusion models[END_REF][START_REF] Ruth | Models of collective cell motion for cell populations with different aspect ratio: diffusion, proliferation and travelling waves[END_REF]. The first model for the diffusion of extended particles with exclusion interactions on a one-dimensional lattice is described in Ref. [START_REF] Schnherr | Exclusion process for particles of arbitrary extension: hydrodynamic limit and algebraic properties[END_REF] for a general process involving symmetric, as well as asymmetric, hopping dynamics of the rods (the theory is named L-ASEP). Referring to rods of length L, where L has to be interpreted as the aspect ratio between the dimensions of the elongated agent, the authors derive a mean-field equation for the one-dimensional exclusion process. In the absence of a field and for equivalent sites the equation for the density of particles reads:

∂ρ ∂t = D 0 ∂ ∂x 2 ρ 1 -(L -1)ρ . (2.87)
The previous equation shows that, at least in one dimension, the extendedcrowding procedure yields a modified diffusion term even in the absence of external fields or spatial inhomogeneities. By defining a density-dependent diffusion coefficient 2 , which has to be regarded as a collective diffusivity, eq. ( 2.87) can be reformulated as a nonlinear transport equation:

D(ρ) = D 0 /[(1 -(L -1)ρ]
∂ρ ∂t = ∂ ∂x D(ρ) ∂ρ ∂x . (2.88)
The nonlinear diffusion equation (2.87) has been derived through an ingenious but complicate change of variables based on a quantitative mapping between the L-ASEP and the zero-range process [START_REF] Schnherr | Exclusion process for particles of arbitrary extension: hydrodynamic limit and algebraic properties[END_REF], however, it turns out that it can be regarded as the local-density approximation (LDA) of a simple general property of one-dimensional exclusion processes. As pointed out in 1967 by Lebowitz and Percus [START_REF] Lebowitz | Kinetic equations and density expansions: Exactly solvable one-dimensional system[END_REF] concerning bulk properties: "For many purposes, however, adding a finite diameter does not introduce any new complications; it merely requires the replacement in certain expressions of the actual volume per particle ρ -1 by the reduced volume ρ -1 -σ, i.e. ρ → ρ/(1 -σρ)."

In that sense, the quantity ρ/[1 -(L -1)ρ] is recognized as an effective density in Ref. [START_REF] Ferreira | Anomalous tag diffusion in the asymmetric exclusion model with particles of arbitrary sizes[END_REF] while computing the velocity of finite sized particles which occupy L units of lattice spacing in one dimension. By performing the substitution

ρ(x, t) → ρ(x, t) [1 -(L -1)ρ(x, t)]
in the Ficks law, one recovers (2.87). Point-like crowding in the mean field approximation corresponds to systems of fully penetrable spheres, while extended-size crowding yields a transport equation suitable for systems of totally impenetrable (hard) spheres. We stress that the case of point-like crowding is recovered for agents of aspect ratio L = 1. The discussion above which leads to equation (2.87) applies to one dimensional systems. Starting from this setting, one can raise the question whether similar arguments might be employed to obtain a modified nonlinear equation accounting for excluded volume effects in the diffusion of hard spheres in two and three dimensions. Unfortunately, the strategy used to recover equation (2.87) in [START_REF] Schnherr | Exclusion process for particles of arbitrary extension: hydrodynamic limit and algebraic properties[END_REF] can not be employed to provide a description of the extended crowding in higher dimensions. Several other models have been proposed for two and three dimensional domains [START_REF] Simpson | Models of collective cell spreading with variable cell aspect ratio: A motivation for degenerate diffusion models[END_REF][START_REF] Penington | Collective motion of dimers[END_REF][START_REF] Bruna | Excluded-volume effects in the diffusion of hard spheres[END_REF], starting from stochastic processes enforced with different microscopic rules. Depending on the shape of the agents (hard rods or hard spheres), on the prescribed hopping rules, on the allowed mechanisms for changing the orientation of the agents (rotation or reptation), and on the mean field assumptions made to recover the macroscopic picture, different equations are derived for the density of extended particles in higher dimensions.

In Ref. [START_REF] Penington | Collective motion of dimers[END_REF] the diffusion of dimers (hard rods with L = 2) is studied from microscopic rules in one and two dimensions. In the one dimensional lattice, for the simplest case of a single population without overlapping processes, the authors propose a discrete model for the probability of a site i being occupied by the right side of an agent after n time steps. A mean-field equation is then obtained by restricting the correlation between the occupancy probabilities to the nearest neighboring sites. The continuum limit for the average right side occupancy R in 1D reads:

∂R ∂t = D 0 ∂ ∂x (1 + 2(L -1)R) ∂R ∂x . (2.89)
The inconsistency between equation (2.87) and (2.89) is only apparent. In fact, due to the approximation on the probabilities of occupancy (the probability of a site being occupied is assumed to depend only on the occupancies of the next (L -1) neighboring sites), eq. (2.89) is more accurate at low densities. At once, the collective diffusivity in (2.88) at low densities becomes D(ρ) ≃ D 0 (1+2(L-1)ρ). For a two-dimensional domain the description is modified by allowing the agents (only dimers are considered here) to be oriented along their longer or shorter dimension, and to change configuration through reptation or traslation. For the simplest nonoverlapping case, after constructing a covering lattice and mapping the problem onto a monomer problem for the motion of the center of the particles, a mean field hypothesis on the correlation of occupancy probabilities for neighboring sites is enforced. The resulting PDE for the total agent occupancy T (both orientations are considered) is:

∂T ∂t = D 0 ∇ • [(1 + 2T )∇T ]. (2.90)
Several extension of this model are considered, for example by accounting for the presence of k different species of dimers, the authors obtain the following equations, respectively in 1D and 2D :

∂R k ∂t = D 0,k ∂ ∂x (1 -R) ∂R k ∂x + 3R k ∂R ∂x , (2.91
)

∂T k ∂t = D 0,k ∇ • (1 - 3 2 T )∇T k + 7 2 T k ∇T (2.92)
where R k is the right side occupancy of species k, R is the right side occupancy of all species, T k is the occupancy of species k and T is the total occupancy of all species.

It is interesting to note that, unlike the single species case (which can be retrieved by summing the equations on the k index), the equation in 1D and 2D for multiple species are different. In Ref. [START_REF] Bruna | Excluded-volume effects in the diffusion of hard spheres[END_REF] the excluded volume effects are analized with regard to the diffusion of finite-size hard speres in a box in two or three dimension. Starting from stochastic

d = 1 d = 2 d = 3 Fully penetrable spheres 1 1 1 Totally impenetrable spheres 1 1 -φ a 0 = 1 (1-φ) 2 b 0 = 1+φ+φ 2 -φ 3 (1-φ) 3
differential equations for the position of the centers of all the N particles, the authors obtain the following equation for the packing fraction φ of the spheres in the limit of large N :

∂φ ∂t (x, t) = ∇ • [(1 + 4(d -1)φ)∇φ] , (2.93)
where d is the space dimension. In all the aforementioned models, for one, two and three dimension, the description of the transport of finite-size impenetrable particles leads to the same structure of the nonlinear diffusion equation. This motivates us to extend the aforementioned arguments by Lebowitz and Percus to higher dimensions. We can identify the substitution made by Lebowitz and Percus as a mapping between certain statistical properties characterizing systems of fully penetrable spheres and totally impenetrable (hard) spheres, with reference to standard definitions of micro-structural descriptors in d dimensions [START_REF] Torquato | Random Heterogenous Materials[END_REF]. More precisely, let us introduce the so-called conditional pair distribution function (CPDF) G p (r). Let r denote the distance from the center of some reference particle in a system with bulk density ρ 0 . Then, by definition ρ 0 s d (r)G p (r) dr equals the average number of particles in the shell of infinitesimal volume s d (r) dr around the central particle, given that the volume v d (r) of the d-sphere of radius r is empty of other particle centers. Here

s d (r) = dv d (r) dr = 2π d/2 r d-1 Γ(d/2) . (2.94)
The CPDF for systems of fully penetrable and totally impenetrable systems are reported in Table (2.1) in the approximation of r ≫ σ which amounts to omit in the definition of the CPDF additional terms which scale as negative powers of the distance. In practical terms, restricting to large values of r amounts to imposing an effective degree of coarse graining, and so neglecting the fine detailed modulations which occur at small spatial scales. More precisely, this choice is expected to be more accurate at low densities, where neglecting the density fluctuations in the first coordinate in a shell does not introduce large errors. Form the analysis of Table 2.1 one can readily recognize the Lebowitz and Percus substitution as a mapping between the FPS and TIS CPDFs in one dimension. This analogy suggests a highly speculative procedure to generalize the L-ASEP to describe excluded-volume effects in more than one dimension in the mean-field approximation in a homogeneous medium and zero field. For a spherically symmetric problem we posit

∂ρ ∂t = D 0 ∇ 2 [ρ Ḡp (ρ)] (2.95)
where it is understood that, according to the local density approximation, one should replace φ with v σ ρ(r, t), where v σ ≡ v D (σ/2) is the volume of one hard d-dimensional sphere, so that Ḡp depends on r implicitly through ρ(r). This is a highly speculative proposal, which however returns the correct description in the limit of low density: operating in this setting, the two-dimensional and three-dimensional non-linear equations for hard spheres (2.93) are recovered from (2.95).

In order to provide another a posteriori justification for eq. (2.95), it is instructive to consider how the classical problem of diffusion to an absorbing sphere is modified in a non-ideal fluid. Let us imagine a fixed sink of radius R s that absorbs hard spheres of radius σ/2 and bulk density ρ 0 . The rate k measuring the number of particles absorbed by the sink per unit time equals the total flux into the sink. For ordinary Fickian diffusion, one has the classical result k .

= k S = 4πD 0 (R s + σ/2)ρ 0 , known as the Smoluchowski rate [START_REF] Smoluchowski | Drei vortrage ubër diffusion brownsche molekular bewegung und koagulation von kolloidteichen[END_REF]. This result is indeed the prediction of a twobody problem, i.e. it amounts to considering the absorption of non-interacting, or equivalently fully penetrable, spheres. Thus, it describes the problem in the infinite dilution limit. Eq. (2.95) can now be employed to repeat the same exercise for hard spheres at finite densities, that is, the Smoluchowski problem with excluded-volume interactions accounted for. One should then solve the following boundary-value problem

∇ 2 [ρ Ḡp (v σ ρ)] = 0 (2.96) ρ(r = R s + σ/2) = 0 (2.97) lim r→∞ ρ(r) = ρ 0 . (2.98)
The rate can be computed readily without really solving the (modified) Laplace equation. From (2.96), we have directly

∂ ∂r ρ(r) Ḡp (φ(r)) = k 4πD 0 r 2 (2.99)
where we have defined φ(r) = v σ ρ(r), so that φ 0 = v σ ρ 0 denotes the bulk packing fraction of the hard spheres. Integrating eq. (2.99) between R s +σ/2 and infinity and taking into account the boundary conditions (2.97) and (2.98), it is straightforward to obtain

k k S = Ḡp (φ 0 ) (2.100)
where Ḡp (φ 0 ) = lim r→∞ Ḡp (φ(r)). In three dimensions, one thus has k/k S = b 0 (φ 0 ), where one can recognize b 0 (φ) as the compressibility Z(φ) of the hard sphere fluid in the Carnahan-Starling approximation [START_REF] Carnahan | Equation of state for nonattracting rigid spheres[END_REF]. We see that eq. ( 2.95) allows one to recover the known result k/k S = Z(φ), obtained in two different ways, by assuming a density-dependent mobility in the diffusion equation [START_REF] Dorsaz | Diffusion-limited reactions in crowded environments[END_REF] and from a transport equation derived in the local-density approximation [START_REF] Piazza | Diffusionlimited reactions in crowded environments: a local density approximation[END_REF].

In all the models derived thus far the diffusion is treated as a non-ideal process because of the high concentration of molecules. In the following part of this thesis, we will abandon the microscopic level to examine a different kind of complexity. We will study the interplay between reaction and diffusion, and the non-ideal conditions will be determined by the complex geometry where the processes are analyzed.

Chapter 3

Theory of diffusion influenced reactions in complex geometries

In the second part of this thesis, we will move on to analyze chemical reactions occurring between different species in complex environments. Aggregation of proteins in composite structures of protein-protein complexes and reactions occurring between microscopic components such as substrates and enzymes are central activities underlying the function of living cells. These processes can be classified as diffusion-influenced reactions (DI), because they involve the diffusion of the reagents and the subsequent chemical fixation of the encounter product. Virtually all biochemical processes in living media can be counted among them, together with those occurring in an ever-growing number of emerging nano-technologies. The general goal of our study is to develop an analytical strategy to compute the reaction rate constant of such processes, while keeping track of the complexity of the environment. In many biological and industrial processes the diffusion of reactants takes place in a confined and crowded environment: cells for example are occupied for over 30% of their volume by membrane-delimited organelles and different sorts of cytoskeletal structures, as shown in Fig. 3.1.

In general, static or moving objects placed in the domain where the reaction occurs can modify the overall behavior of the diffusing agents and the efficiency of the binding process. Inside a cell, for example, a diffusing ligand encounters a large number of sinks, traps and obstacles that can slow down or even arrest its motion towards the nucleus or the target receptor in the cytoplasm. The extent of the alteration depends on the characteristic of the neighboring elements: number, size, position, shape, and their nature either reactive or non-reactive. To simplify the analysis of the reaction kinetics, these effects are often neglected for systems which are enough diluted to allow one to overlook the influence of other reactants and other space-filling species. In this limit the DI process can be represented as a two-body reaction and tackled through well known analytical and numerical schemes. On the other hand, in the majority of biological and technological media, the concentration of receptors, obstacles and reactive boundaries is sufficiently high to require an accurate analysis of the many-body interactions.

A general way to describe the formation of a product from a bimolecular reaction between two microscopic components A and B in solution is:

[A] + [B] k D ---⇀ ↽ --- k -D [AB] * k * -→ [P ] (3.1)
where, following the standard notation, the symbol [•] denotes the concentration. For example, this might represent the catalysis of a substrate-product conversion or the formation of a stable protein-protein complex. The first step leads to the formation of the encounter complex [AB] * . This process is driven by diffusion and characterized by a rate constant k D . Once the encounter complex is formed, it can be either chemically fixed to give rise to a stable complex (with a finite probability expressed by k * , or else the two molecules can diffuse away (with rate constant k -D ).

In the quasi-equilibrium approximation, which amounts to consider a stationary intermediate product concentration (or equivalently, to be in a situation where the concentration of the substrate molecules B that are turned into product P is much larger than the concentration of A), the reaction scheme (3.1) can be approximated by the pseudo first-order reaction

[A] + [B] k -→ [P ] (3.2) with reaction rate constant k = k D k * (k * + k -D )
. k is the rate of formation of products, i.e. the number of B molecules turned into product P per unit time (it has the physical dimension of [t -1 ]).

For the case of chemical reactions connoted by a high chemical affinity, such as antigen-antibody reactions (see section 5), one has k * ≫ k -D and the limiting step for the process is given by diffusion, because the agents react immediately upon getting in contact. The overall rate constant in this case is equal to the encounter rate k ≃ k D . The reactions belonging to this class are named diffusion-controlled or diffusion-limited.

In 1917 M. Von Smoluchowski showed how to compute the rate constant k D that controls the diffusion-limited formation of the encounter complex between a ligand (B) and a receptor (A) in a simple approximation of high dilution. The idea is that, under the hypotheses that:

• both ligand and receptor are represented as diffusing spherical bodies

• A molecules diffuse much more slowly than B molecules (they are of much larger size),

• both species are highly diluted

• the bulk concentration of A molecules is much smaller than the bulk concentration of B molecules, it is possible to consider a single receptor at rest in the domain where the ligands are diffusing, and the chemical reaction can be mapped onto a two-body boundary problem [START_REF] Sanchez-Osorio | Foundations for modeling the dynamics of gene regulatory networks: a multilevel-perspective review[END_REF]. More precisely, the receptor is represented as a spherical sink centered in the reference frame, with radius equal to the sum of the linear dimensions of the two molecules, R = R A + R B , while the ligands are point-like agents diffusing in the open domain with the effective diffusivity D = D A + D B . When a molecule of the species B reaches the boundary it is immediately absorbed: this physically amounts to say that the two agents react immediately when they get in contact. Here the hypothesis of high dilution allows one to neglect the encounter probability of two A and the effects of crowding for the ligands, in the sense that two B agents never interact with each other. With those assumptions, in the long time limit (we are assuming t >> t D , being t D the characteristic relaxation time for diffusion), the concentration of the diffusing population obeys the Laplace equation, enforced with absorbing boundary conditions at the encounter surface, namely:

     ∇ 2 c = 0 c| |r|=R = 0 lim |r|→∞ c = c 0 . , (3.3)
where c 0 is the bulk concentration of ligands and the boundary |r| = R is a spherical sink of effective radius R A + R B (the ligand-receptor encounter distance). This description applies for example to a ligand that is diffusing towards the receptorcovered surface of a cell (see figure 3.4 and the discussion in the conclusive chapter).

The solution of the stationary problem can be worked out explicitly as:

c(r) = c 0 1 - R r . (3.4)
The overall rate of the reaction (the parameter k defined in 3.2) corresponds mathematically to the flux across the sink surface (number of B molecules crossing the surface per unit time):

k = - ∂Ω J • n dσ = D ∂Ω ∂c ∂r dσ = 4πDRc 0 := k S
where n is the normal vector to the surface of integration and ∂Ω is the boundary of the sink. The specific value of the rate constant relative to the Smoluchowski model is customarily denoted as k S1 . The perfectly absorbing boundary condition at the separation surface of the two reactants characterizes the process as diffusion limited, because the reaction happens with probability 1 when the agents are at the encounter distance.

Our goal is to extend this simple setting by relaxing the hypotheses of extreme dilution and open domain, and accounting for partially reactive surfaces and more complex shapes of receptors. Before moving on to consider the many-body problem, it is useful to spend a few words to see how to introduce the partial reactivity in the mathematical description of the simple two-body setting (3.3). The idea, first proposed by Collins and Kimball in 1949 in [START_REF] Frank | Diffusion-controlled reaction rates[END_REF], is that the chemical fixation rate k * can be thought of yielding a finite absorption probability for ligand molecules: the encounter complex only gives rise to the chemical transformation with a certain

c(r) = c 0 1 - R r k * k * + 4πDR
and yields the modified rate constant:

k = k S k * k * + 4πDR
.

By taking k * → ∞ in the previous expression, we recover the Smoluchowski rate constant k S , while in the limit k * = 0, the rate constant vanishes and the BCs become perfectly reflecting, ∂c/∂r = 0. In the chemical physics community, chemical reactions which are neither rate-limited by diffusion nor by reaction, that is with intermediate values of k * , are termed generically diffusion-influenced [START_REF] Szabo | Theory of diffusion-influenced fluorescence quenching[END_REF]. Now that we have shown how to account for the partial reactivity of a given boundary, we can move on to describe the multi-body problem, which will provide a more suitable instrument for modeling diffusion-reaction processes in complex environments. In fact, in many realistic situations in chemical and biochemical kinetics, a single ligand (B) has to diffuse among many competing reactive particles (A). In addition, it might be forced to find its target within a specific confining geometry, which in principle can be modeled through a collection of reflecting boundaries.

Let us consider a confined spherical domain, containing an arbitrary number of spherical boundaries, each one characterized by a specific size and a specific reactivity. These structures are meant to represent all the reactive or inert bodies with which a diffusing molecule gets in contact while moving in a complex medium. For example, in describing intracellular diffusion, cellular organelles and other cytoskeletal structures which do not react with the specific ligands, can be represented as fully reflecting spheres or assemblies of spheres to make up other non-spherical shapes. This choice amounts to consider k * = 0 in the BC, while the reactive targets are associated with values of k * > 0. In this way, a diffusing molecule which gets in contact with one of the hard spheres can either be absorbed if the body is reactive, or just have its motion obstructed due to the steric effects of a nonreactive body.

To match this description, the Laplace equation for the density of ligands is enforced with (N +1) radiative boundary conditions, where Ω 1 , . . . , Ω N is the collection of reactive boundaries included in the domain Ω 0 :

             ▽ 2 c = 0 4πDR 2 α ∂c ∂r α -k * α c ∂Ωα = 0 ∀α = 1...N ∂c ∂r 0 -h 0 (c 0 -c) ∂Ω 0 = 0, . (3.6)
Here c 0 stands for the constant bulk density of ligands, k * α , R α are the chemical and geometrical parameters which characterize the selected boundary, and h 0 specifies the reactivity of the inner surface of the confining sphere Ω 0 (we will discuss its interpretation in subsection 3.1.2). The solution of this boundary problem can not be worked out as easily as in the previous cases, since the global spherical symmetry is lost by taking into account multi-body boundary conditions. We will see in the next section that the confinement of the domain and the presence of multiple spheres require the use of more elaborate mathematical tools.

Analytical solution

The problem stated in the previous section can be solved both in the unbounded domain outside the given set of spherical boundaries, or within a bounded spherical domain that contains all the reactive boundaries. For the moment, we will restrict our analysis to the domain:

Ω = Ω 0 \ N α=1
Ω α as represented in Fig. 3.4. In the following (see chapter 4) we will show how to describe a situation where we consider ∂Ω 0 (the boundary of Ω 0 ) as a permeable barrier between two regions with different properties.

Let us introduce the non dimensional quantities: the density normalized on the bulk value u(x) = c(x)/c 0 and the variables ξ α = r α /R α and ξ 0 = r 0 /R 0 , normalized on the radii of the reactive boundaries. The boundary problem (3.6) becomes:

             ▽ 2 u = 0 ∂u ∂ξ α -h α u ∂Ωα = 0 ∀α = 1...N ∂u ∂ξ 0 -h 0 (1 -u) ∂Ω 0 = 0. (3.7)
We will refer to the parameter h α = k * α /4πDR α , which determines the reactivity of the α-th sphere, as to the intrinsic reactivity of the selected boundary. Similarly, the parameter h 0 is the intrinsic reactivity of the inner surface of the spherical domain Ω 0 , which can be interpreted as carrying some information on the space external to ∂Ω 0 , even if the problem is considered only inside. More specifically, it can be interpreted as the ratio between the diffusion constants within the media lying at either ) is employed to model intracellular diffusion, and ∂Ω 0 is meant to represent the cell membrane this ratio is usually < 1. Boundary problems of the kind (3.7), featuring multi-connected domains, are customarily solved numerically, with finite element methods, employed to avoid the substantial analytical difficulties of the calculations. Even if the numerical strategy proves to be effective, due to its high versatility, it does not provide any analytical insight and it involves extremely high computational costs, even for a limited number of reactive boundaries.

Our aim is thus to provide an analytical solution, which enables one to solve the problem to any desired accuracy for an arbitrary composition of the domain, which means for any choice of number, position and dimension of the reactive boundaries. The idea is to take into account (N + 1) spherical coordinate systems, one for each spherical boundary considered in the domain. Then we can write the solution of the Laplace equation as a sum of linear combinations of spherical harmonics, each one expressed in the spherical system associated with the selected αth boundary. The global solution, which is formally expressed in (N + 1) different coordinate systems, contains one set of coefficients per each harmonic expansion. The global set of coefficients must be determined as usual by imposing the boundary conditions. This poses the biggest challenge, because in the neighborhood of each sphere this operation requires to express the global density as a function of the local coordinates. This is required since each boundary equation is expressed in spherical symmetry with respect to the corresponding local coordinate system.

Analytical solution 55

We start by considering the (N + 1) spherical coordinate systems corresponding to each boundary (see figure (3.1)): r 0 = (r 0 , θ 0 , φ 0 ) and r α = (r α , θ α , φ α ) ∀α = 1, 2, . . . N , where θ and φ identify, respectively, the polar and the azimuthal angles. It is known that the solution of the Laplace equation within a spherical boundary can be written as a linear combination of regular harmonics:

u + 0 = ∞ n=0 n m=-n A mn u + mn (r 0 ) = ∞ n=0 n m=-n A mn ξ n 0 Y mn (r 0 ) = ∞ n=0 n m=-n A mn r 0 R 0 n P m n (cos θ 0 )e imφ 0 , (3.8)
while a basis of irregular harmonics must be used to express the solution within an unbounded domain outside a given spherical boundary (the regular harmonics diverge at infinity):

u - α = ∞ n=0 n m=-n B α mn u - mn (r α ) = ∞ n=0 n m=-n B α mn ξ -n-1 i Y mn (r α ) = ∞ n=0 n m=-n B α mn r α R α -n-1 P m n (cosθ α )e imφα .
(3.9)

Using the superposition principle for the Laplace equation, we can write the solution in Ω as a sum of linear combinations of regular (for solving (3.7) inside Ω 0 ) and irregular harmonics (for expressing the solution outside each Ω α ):

u = u + 0 + N α=1 u - α = ∞ n=0 n m=-n A mn ξ n 0 Y mn (r 0 ) + ∞ n=0 n m=-n B α mn ξ -n-1 i Y mn (r α ).
The values of coefficients of the above expansion are determined by imposing the boundary conditions. In the neighborhood of each boundary we have to express all the bases as a function of the local coordinates. More precisely, to determine the values of A mn and B α mn , in a neighborhood of each ∂Ω α (α = 1, 2, . . . , N ), we have to express u + 0 and u - j , j = i, as a function of the r α coordinates, and similarly, in a neighborhood of ∂Ω 0 , we have to write every u - α as a function of r 0 . To this purpose, we will make use of the addition theorems for spherical harmonics, [START_REF] Morse | Methods of Theoretical Physics, Part II[END_REF].

Let us start by imposing the condition on ∂Ω 0 . Let L 0α be the constant vector connecting the center of ∂Ω 0 to the center of ∂Ω α , so that r α = r 0 -L 0α . The relation |L 0α | < |r 0 | holds for all value of α = 1, 2, . . . , N in a neighborhood of ∂Ω 0 , thus for every u - α we need the re-expansion formula from irregular to irregular harmonics (third panel on fig. 3.6)

r -n-1 α Y mn (r α ) = ∞ l=0 n s=-n (-1) l+s (n + l -m + s)! (n -m)!(l + s)! L l 0α Y sl (-L 0α )r -(n+l)-1 0 Y m-s,n+l (r 0 ) (3.10)
that after the change of indexes:

n + l = q m -s = g (3.11)
Theory of diffusion influenced reactions can be written as:

r -n-1 α Y mn (r α ) = ∞ q=n m+(q-n) g=m-(q-n) (-1) q-n+m-g (q -g)! (n -m)!(q -n + m -g)! L q-n 0α Y m-g,q-n (-L 0α )r -q-1 0 Y gq (r 0 ).
(3.12)

After multiplying by R n+1 α we obtain the non-dimensional form of the re-expansion formula:

ξ -n-1 α Y mn (r α ) = ∞ q=n m+(q-n) g=m-(q-n) (-1) q-m+n-g (q -g)! (n -m)!(l + s)! L 0α R 0 q-n Y m-g,q-n (-L 0α ) R 0 R α -n-1 ξ -q-1 0 Y gq (r 0 ) = ∞ q=n m+(q-n) g=m-(q-n) V α,m,n g,q ξ -q-1 0 Y gq (r 0 ) (3.13)
where

V α,m,n g,q = (-1) q-n+m-g (q -g)! (n -m)!(q -n + m -g)! L 0α R 0 q-n R 0 R α -n-1 Y m-g,q-n (-L 0α ) = (-1) q-m+n-g (q -g)! (n -m)!(q -n + m -g)! η q-n 0α χ n+1 α Y m-g,q-n (-L 0α ). (3.14)
and we have defined ∀i, j = 0, 1, . . . , N :

η i,j := L i,j R 0 with η i,j = η j,i χ i := R i R 0 . (3.15)
Hence the solution of our problem in the vicinity of ∂Ω 0 can be written as

u = ∞ n=0 n m=-n A mn ξ n 0 Y mn (r 0 ) + N α=1 ∞ n=0 n m=-n B α mn ∞ q=n m+(q-n) g=m-(q-n) V α,m,n g,q ξ -q-1 0 Y gq (r 0 ) = ∞ q=0 q g=-q A gq ξ q 0 Y gq (r 0 ) + N α=1 ∞ n=0 n m=-n B α mn ∞ q=n m+(q-n) g=m-(q-n) V α,m,n g,q ξ -q-1 0 Y gq (r 0 ) = ∞ q=0 q g=-q Y gq (r 0 ) A gq ξ q 0 + N α=1 ∞ n=0 n m=-n B α mn V α,m,n g,q ξ -(q+1) 0 1 {q≥n} 1 {m-(q-n)≤g≤m+(q-n)} = ∞ q=0 q g=-q Y gq (r 0 ) A gq ξ q 0 + N α=1 q n=0 n m=-n B α mn V α,m,n g,q ξ -(q+1) 0 1 {g-(q-n)≤m≤g+(q-n)} (3.16)
3.1 Analytical solution 57

while the first derivative of the stationary density field with respect to ξ 0 reads:

∂u ∂ξ 0 = ∞ q=0 q g=-q Y gq (r 0 ) qA gq ξ q-1 0 -(q + 1)ξ -(q+2) 0 N α=1 q n=0 n m=-n B α mn V α,m,n g,q 1 {g-(q-n)≤m≤g+(q-n)} .
(3.17)

In the previous expression we have introduced the characteristic function:

1 A = 1 A (x) = 1 if x ∈ A 0 if x / ∈ A (3.18)
where the subset A is identified by a set of inequalities involving the indices of the spherical harmonics. Having obtained the expression of the density in the coordinate system relative to Ω 0 allows us to impose the boundary condition on Ω 0 , namely:

∞ q=0 q g=-q Y gq (r 0 ) qA gq ξ q-1 0 + (-q -1)ξ -(q+2) 0 N α=1 q n=0 n m=-n B α mn V α,m,n g,q 1 {g-(q-n)≤m≤g+(q-n)} -h 0 A gq ξ q 0 + N α=1 q n=0 n m=-n B α mn V α,m,n g,q ξ -(q+1) 0 1 {g-(q-n)≤m≤g+(q-n)} + h 0 ∂Ω 0 = 0. (3.19)
Since the spherical harmonics are a base on the sphere, equation (3.19) implies:

A gq (q -h 0 ) + (-h 0 -q -1)

N α=1 q n=0 n m=-n B α mn V α,m,n g,q 1 {g-(q-n)≤m≤g+(q-n)}
= -h 0 δ ((g,q),(0,0)) (3.20) which gives for every q = 0, 1, . . . , N and g = -q, -q + 1, . . . , q -1, q:

A gq + -h 0 -q -1 (q -h 0 ) N α=1 q n=0 n m=-n B α mn V α,m,n g,q 1 {g-(q-n)≤m≤g+(q-n)} = 1. (3.21)
Let us now consider the boundary conditions on the N spheres within Ω 0 . In the neighborhood of each ∂Ω α , we must express u + 0 and each u - β with β = α in the r α coordinates. Let us first focus on the re-expansion of u + 0 . Once we have expressed r 0 = r α +L 0α , where L 0α is the constant vector connecting the center of Ω 0 with the center of Ω α , we can use the re-expansion formula from regular to regular harmonics (see first panel on fig. 3.6)

r n 0 Y mn (r 0 ) = n q=0 q g=-q (n + m)! (n -q + m -g)!(q + g)! L n-q 0α Y m-g,n-q (L 0α )r q α Y gq (r α ). (3.22)

Theory of diffusion influenced reactions

After dividing by R n 0 we obtain the non-dimensional relation:

ξ n 0 Y mn (r 0 ) = n q=0 q g=-q n + m q + g L 0α R 0 n-q r α R α q R α R 0 q Y m-g,n-q (L 0α )Y gq (r α ) = n q=0 q g=-q H (α,g,q) m,n ξ q α Y gq (r α ) (3.23)
where we have defined the matrix H:

H (α,g,q) m,n = n + m q + g χ q α η n-q 0α Y m-g,n-q (L 0α ).
In order to express u - β with β = α on ∂Ω α , we consider r β = L βα +r α . Since the relation |L βα | > |r α | holds in a vicinity of ∂Ω α , we must use the irregular-to-regular re-expansion formula (see the second panel in fig. 3.6)

r -n-1 β Y mn (r β ) = ∞ q=0 q g=-q (-1) q+g (n -m + q + g)! (n -m)!(q + g)! L -(n+q)-1 βα Y m-g,n+q (L βα )r q α Y gq (r α ).
(3.24)

that becomes as a function of dimensionless variables:

ξ -n-1 β Y mn (r β ) = ∞ q=0 q g=-q (-1) q+g (n -m + q + g)! (n -m)!(q + g)! L βα R 0 -(n+q)-1 R α R 0 q r α R α q R 0 R β -n-1 Y m-g,n+q (L βα )Y gq (r α ) = ∞ q=0 q g=-q (-1) q+g (n -m + q + g)! (n -m)!(q + g)! η -(n+q)-1 βα χ q α χ n+1 β Y m-g,n+q (L βα )ξ q α Y gq (r α ) = ∞ q=0 q g=-q W (α,β,g,q)
m,n ξ q α Y gq (r α ).

(3.25)

where

W (α,β,g,q) m,n = (-1) q+g (n -m + q + g)! (n -m)!(q + g)! η -(n+q)-1 βα χ q α χ n+1 β Y m-g,n+q (L βα ). (3.26)
Hence in the vicinity of each ∂Ω α we can write the solution as :

u = ∞ n=0 n m=-n B α mn ξ -n-1 α Y mn (r α ) + N β=1,β =α ∞ n=0 n m=-n B β mn ξ -n-1 β Y mn (r β ) + ∞ n=0 n m=-n A mn ξ n 0 Y mn (r 0 ) = ∞ q=0 q g=-q B α gq ξ -q-1 α Y gq (r α ) + N β=1,β =α ∞ n=0 n m=-n B β mn ∞ q=0 q g=-q W (α,β,g,q) m,n ξ q α Y gq (r α )+ ∞ n=0 n m=-n A mn n q=0 q g=-q H (α,g,q) m,n ξ q α Y gq (r α )1 q≤n = ∞ q=0 q g=-q Y gq (r α ) B α gq ξ -q-1 α + N β=1,β =α ∞ n=0 n m=-n B β mn W (α,β,g,q) m,n ξ q α + ∞ n=q n m=-n
A mn H (α,g,q) m,n ξ q α .

(

The derivative with respect to ξ α is:

∂u ∂ξ α = ∞ q=0 q g=-q Y gq (r α ) (-q -1)B α gq ξ -q-2 α + qξ q-1 α N β=1,β =α ∞ n=0 n m=-n B β mn W (α,β,g,q) m,n + ∞ n=q n m=-n
A mn H (α,g,q)

m,n .

(3.28)

Imposing the boundary condition yields :

∞ q=0 q g=-q Y gq (r α ) (-q -1 -h α ξ α )B α gq ξ -q-2 α + (qξ q-1 α -h α ξ α ) N β=1,β =α ∞ n=0 n m=-n B β mn W (α,β,g,q) m,n + ∞ n=q n m=-n A mn H (α,g,q) m,n ∂Ωα = 0, (3.29)
which gives for every q = 0, 1, . . . , N and g = -q, -q + 1, . . . , q -1, q:

-(q + 1 + h α )B α gq + (q -h α ) ∞ n=0 n m=-n (A mn H (α,g,q) m,n 1 q≤n + N β=1,β =α B β mn W (α,β,g,q) m,n ) = 0.
(3.30) We have thus obtained the second equation:

-B α gq + (q -h α ) (h α + q + 1) ∞ n=0 n m=-n A mn H (α,g,q) m,n 1 q≤n + N β=1,β =α B β mn W (α,β,g,q)
m,n = 0.

(3.31) To summarize, by using the addition theorems together with the uniqueness of the spherical harmonics expansions, we managed to map the boundary problem (3.7) to the solution of an infinite-dimensional linear system that has to be solved for the coefficients of the spherical harmonics expansion. More precisely, we have for all the multipoles q = 0, 1, . . . , N and g = -q, -q + 1, . . . , q -1, q and for each reactive boundary Ω α with α = 1, 2, . . . , N :

             -B α gq + (q -h α ) (h α + q + 1) ∞ n=0 n m=-n A mn H (α,g,q) m,n 1 q≤n + N β=1,β =α B β mn W (α,β,g,q)
m,n = 0

A gq + -h 0 -q -1 (q -h 0 ) N α=1 q n=0 n m=-n B α mn V α,m,n g,q 1
{g-(q-n)≤m≤g+(q-n)} = δ (g,q)=(0,0) .

(3.32) where

V α,m,n g,q = (-1) q-n+m-g (q -g)! (n -m)!(q -n + m -g)! η q-n 0α χ n+1 α Y m-g,q-n (-L 0α ) (3.33) H (α,g,q) m,n = n + m q + g χ q α η n-q 0α Y m-g,n-q (L 0α ) (3.34) W (α,β,g,q) m,n = (-1) q+g (n -m + q + g)! (n -m)!(q + g)! η -(n+q)-1 βα χ q α χ n+1 β Y m-g,n+q (L βα ), (3.35)
where

χ i := R i R 0 and η i,j := L ij R 0 with η ij = η ji .
The matrix representation of system (3.32) reads:

                           1 V 1 V 2 . . . V N H 1 -1 W 1,2 . . . W 1,N H 2 W 2,1 -1 . . . W 2,N . . . . . . . . . . . . . . . H N W N,1 W N,2 . . . -1                            ×                            A 00 . . . A N M N M B 1 00 . . . B 1 N M N M . . . B N 00 . . . B N N M N M                            =                            1 . . . 0 0 . . . 0 . . . 0 . . . 0                           
where the physical, chemical and geometrical parameters (radii of the sinks, distances between the centers, reactivities of the boundaries) are hidden inside the matrices W, V, H.

In order to solve the infinite-dimensional system (3.32) one needs to truncate the expansion, by including a finite number of multipoles, N t , which ensures that a desired accuracy is attained. The resulting truncated system for N internal boundaries comprises (N + 1)(N t + 1) 2 equations.

The rate of a reaction between the ligand and a selected boundary Ω α can be computed easily, by calculating the incoming flux as

k α = - ∂Ωα Jds where J = -D ∂ρ ∂r 1 = D ρ 0 R 1 ∂u ∂ξ 1 .
Substituting equation (3.28) in the definition of k α , and using the properties of

Theory of diffusion influenced reactions

Legendre polynomials , we get:

k α k Sα = 1 2 ∞ q=0 q g=-q (-q -1)B α gq + q N β=1,β =α ∞ n=0 n m=-n B β mn W (α,β,g,q) m,n + ∞ n=0 n m=-n A mn H (α,g,q) m,n 1 q≤n 2π 0 e igφα dφ α 1 -1 P gq (µ α )dµ α = 1 2 ∞ q=0 (-q -1)B α 0q + q N β=1,β =α ∞ n=0 n m=-n B β mn W (α,β,0,q) m,n + ∞ n=0 n m=-n
A mn H (α,0,q) m,n 1 q≤n

1 -1 P 0q (µ α )dµ α = -B α 00 (3.36)
where k Sα = 4πDR α c 0 is the Smoluchowski rate computed for an isolated sphere of the same radius in the unbounded domain. The procedure that has led us to the result is quite cumbersome, but we stress that once the initial effort is made to derive the structure of the linear system, we can directly use it to compute the reaction rate to an arbitrary distribution of boundaries. More precisely, given an arrangement of spherical boundaries, it suffices to substitute the structural and chemical parameters that characterize the configuration into the matrices V, W, H and then choose a tolerance for the solution.

The analytical method that we have developed is an useful tool. On the one hand it can be employed to compute the reaction rate constant for a broad choice of reactive landscapes, made by assembling spherical boundaries of selected size at given locations in space and endowed with arbitrary surface reactivity, with much less computational effort than required by the usual numerical methods. On the other, the linear system (3.32) allows one to derive simple approximate analytical formulas that can be used to investigate naturally occurring reactive geometries and to assist in the design of artificial technological setups such as core-shell or yolkshell nanoreactors. In ref. [START_REF] Gordelyi | Transient heat conduction in a medium with multiple spherical cavities[END_REF] a similar method based on expansion formulas is used for solving the time dependent problem associated with eq. (3.7). The authors explicitly compare the efficiency of the analytical calculation with numerical methods based on finite elements, showing how the analytic approach reduces the computation time from the orders of hours to the order of minutes. We note that our method is utterly general, as it can be easily extended to accommodate for reactive environments realized with more complex, non-spherical boundaries. The only requirement is that one of the coordinate systems in which Laplace's equation is separable be used [START_REF] Mccord | Methods of theoretical physics[END_REF], and that addition theorems exist for the corresponding elementary solutions [START_REF] Caola | Solid harmonics and their addition theorems[END_REF].

Rules for selecting the appropriate Addition Theorem

The addition theorems for spherical harmonics allow one to express a combination of spherical harmonics, written in multiple coordinate systems, as a function of a chosen one. They were of fundamental importance for the analytical derivation of the system (3.32), because they allowed us to impose the boundary condition on each sphere. Depending on the type of spherical harmonic that one wants to reexpand (regular or irregular) and on the geometrical parameters characterizing the

r n β Y mn (r β ) = n q=0 q g=-q (n + m)! (n -q + m -g)!(q + g)! L n-q βα Y m-g,n-q (L βα )r q α Y gq (r α ). (3.37)
If one has to re-expand an irregular harmonic u -(r β ), two cases are possible, depending on the ratio between the distance L β α between the centers of the old and new reference systems, and the norm of the the vector r α expressing the position of P in the new system S α . More precisely, if |r α | < |L βα |, then we have to write the irregular harmonic as a function of the regular harmonics centered in S α : To summarize, one can use the following scheme to change variable from system S β to S α :

r -n-1 β Y mn (r β ) = ∞ q=0 q g=-q (-1) q+g (n -m + q + g)! (n -m)!(q + g)! L -(n+q)-1 βα Y m-g,n+q (L βα )r q α Y gq (r α ). ( 3
• u + (r β ) = f (u + (r α )) • u -(r β ) = f (u + (r α )) if |r α | < |L βα | f (u -(r α )) if |r α | > |L βα | .

An interpretation for the boundary condition on the internal boundary Ω 0

In this section, we provide a plausible physical interpretation for the boundary conditions imposed on the internal surface of the spherical cavity Ω 0 considered in the problem depicted in Fig. 3.8. Let us consider the following stationary problem defined in

R 3 \ Ω 1 in spherical coordinates          ∇ • D(r)∇ρ = 0 ρ| r=R 1 = 0 lim r→∞ = ρ B (3.40)

Axially symmetric problems 65

where the diffusion coefficient is a step function

D(r) = D in if R 1 < r < R 2 D out if r > R 2 .
The solution can be found as

ρ(r) = ρ B r r=R 1 dx x 2 D(x) ∞ R 1 dx x 2 D(x) -1 , (3.41)
and the rate of the reaction k becomes

k = 4πR 2 1 D in ∂ρ ∂r r=R 1 = 4πρ B ∞ R 1 dx x 2 D(x) -1 = k S (1 -ǫ) + χǫ (3.42) where χ = D in /D out and ǫ = R 1 /R 2 and k S = 4πD in R 1 ρ B . Given ǫ = R 1 /R 2 < 1,
and the fact that for many applications the ratio χ = D in /D out < 1, the rate of the reaction is greater than the Smoluchowski rate to a sink embedded in an infinite domain characterized by a diffusivity D in :

k = k S 1 -ǫ(1 -χ) > k S . (3.43)
We want to understand which are the hypotheses that make the system (3.40) equivalent to the following boundary value problem considered in the bounded domain

Ω 2 \ Ω 1          ∇ 2 ρ = 0 ρ| r=R 1 = 0 4πD in R 2 2 ∂ρ ∂r -k * (ρ B -ρ) r=R - 2 = 0, (3.44)
where k * is the intrinsic rate constant. After defining h = k * /4πD in R 2 , it is straightforward to show that the function

ρ(r) = ρ B h h(1 -ǫ) -ǫ 1 - R 1 r (3.45)
solves the problem (3.44) in Ω 2 \ Ω 1 , and yields the reaction rate

k = h h(1 -ǫ) + ǫ . (3.46)
By comparing the reaction rates (3.42) and (3.46) obtained for the two systems, we see that the two problems are equivalent provided that

h = 1 χ = D out D in . ( 3

Axially symmetric problems

It is instructive to consider axially symmetric problems as a start to illustrate our methods, as the calculations are substantially less involved. Let us consider an axially symmetric domain, with the centers of the spheres Ω 0 , Ω k lying on the same line, as depicted in figure (3.9). The origin of our reference system is located at the center of Ω 0 with the z-axis passing through the centers of the other spheres. This choice allows us to take into account only the variables r i and θ i because the solution is no longer dependent on the azimuthal angles φ i . The stationary concentration of diffusing molecules for a configuration involving N internal spheres centered on the z-axis is simply obtained by setting g = m = 0 in the system (3.32), since the solution is independent of φ. For this reason, the expansions involve the simpler expressions

u + 0 = ∞ n=0 A n u + n (r 0 ) = ∞ n=0 A n ξ n 0 P n (cosθ 0 ) u - α = ∞ n=0 B n u - n (r - α ) = ∞ n=0 B n ξ -n-1 α P n (cosθ α ).
(3.48)

Axially symmetric problems 67

In this case, the system (3.32) reduces to

           -B α q + (q -h α ) (h α + q + 1) ∞ n=0 A n H (α,0,q) 0,n 1 q≤n + N β=1,β =α B β n W (α,β,0,q) 0,n = 0 A q + -h 0 -q -1 (q -h 0 ) N α=1 q n=0 B α n V α,0,n 0,q = δ (q,0) , (3.49)
with

V α,0,n 0,q = (-1) q-n q n L 0α R 0 q-n R α R 0 n+1 P q-n (cos θ L α0 ), (3.50)
H (α,0,q) 0,n = n q R α R 0 q L 0α R 0 n-q P n-q (cos θ L 0α ), (3.51)
W (α,β,0,q) 0,n = (-1) q n + q q L βα R 0 -n-q-1 R α R 0 q R β R 0 n+1 P n+q (cos θ L βα ). (3.52)
The structure of the matrices H, V , W is much simpler here because the Legendre polynomials are evaluated in 1 or -1 (we recall that P n (-1) = (-1) n , P n (1) = 1) because θ L βα = ±π depending on the relative position of the spheres. Moreover the dimension of the linear system (3.49) is substantially reduced: fixing a number N M of multipoles leads to a system of (N + 1)N M equations.

The concept of diffusive interactions

When a ligand can react with two neighboring boundaries, the overall reaction rate is smaller than the sum of the rates associated with two isolated boundaries. This effect amounts to the existence of an effective interaction between two neighboring reactive surfaces. The described interaction goes by the name of Diffusive Interaction (DI). To illustrate this, let us consider two perfectly absorbing sinks whose centers are separated by a distance L. In 1976 Deutch, Felderhof and Saxton first introduced the concept of DI in the context of this problem [START_REF] Deutch | Competitive effects in diffusioncontrolled reactions[END_REF], whose exact solution was first reported by Samson and Deutch in [START_REF] Samson | Exact solution for the diffusion controlled rate into a pair of reacting sinks[END_REF]. In the context of our theory, an approximate solution can be obtained easily within the monopole approximation (MOA) which amount to keeping but the q = 0 and n = 0 terms (N t = 1) in the system (3.49). By doing this in the unbounded domain (R 0 → ∞) , it is straightforward to show that the rate per sink in the two-sinks problem reads [START_REF] Deutch | Competitive effects in diffusioncontrolled reactions[END_REF]

k k S = 1 1 + R/L . (3.53)
As mentioned before, the overall reaction rate is less than 2k S , i.e. less than the sum of the fluxes into two isolated sinks. The two spheres compete for the flux of ligand molecules and shield part of the flux from each other. Furthermore, one immediately notices that the diffusive interaction is a long range effect, akin to the Coulomb interaction since it persists within any finite distance. This result was extended by S. Traytak in ref. [START_REF] Traytak | The diffusive interaction in diffusion-limited reactions: the steadystate case[END_REF] to the case of two different sinks of radii R α and R β :

k α k Sα = 1 -ǫ β 1 -ǫ α ǫ β where ǫ α = R α L , ǫ β = R β L . ( 3
k α k Sα = 1 -R/L 1 -c R/L 2 . (3.55)
If we analyze the behavior of k α /k Sα as a function of the distance L between the centers, we find a non monotonic trend for selected values of the ratio c. More specifically, a minimum for the rate constant exists at

L = cR + R √ c 2 -c for values of c > (1 + √ 5)/2.
This means that while the rate of the smaller sink increases monotonically with the distance L from the other particle, the same value for the larger sphere is first seen to decrease for small L (see Fig. 3.10). This is the result of the competition between two effects. When the small particle lies very close to the surface of the large sink Ω α , the latter behaves as an effective isolated sink of size R ef f < R α , absorbing a flux k ef f = 4πDR ef f . Upon increasing the distance L, the total flux to the small sink will increase (its active surfaces get larger). The effect of this on the flux to Ω α will depend on the size of the particle Ω β . If R β is small enough, R ef f is not much smaller than R α , so that its flux at the contact distance is not much smaller than k α (L → ∞) = k Sα . Under these conditions, the flux into the large sink starts decreasing, as the other screening particle effectively steal more and more flux from it. However, upon increasing L past a critical distance, the small particle can no longer catch enough ligand flux, so that the flux to Ω α starts increasing, as it should, towards k Sα . Normalized rate constant The position of the minimum, that is absent for c = 1, is shifted to the right as c is increased.

δ / R c = 1 c = 2 c = 3 c = 4
To verify that the non-monotonic growth is a physical effect and not just an anomaly introduced by the monopole truncation, we computed the exact reaction rate by solving the system (3.49) and we compared it with the monopole approximation. As it is clear from Fig (3.11), the MOA provides a good interpolation of the exact results, and the accuracy increases with the the ratio between the sinks' sizes. The position of the minimum is also well approximated by eq. (3.55).

As a second simple illustration of our theoretical procedure, let us consider the rate into a sink confined within a spherical cavity bonded by a permeable membrane. This problem has been considered in ref. [START_REF] Chen | Diffusion into a nanoparticle with first-order surface reaction confined within a sphere[END_REF] where a solution is worked out by posing the calculation in bispherical coordinates. This analysis is motivated by the fact that some bimolecular reactions are enhanced by confining the diffusion of the reagents to a restricted domain (see for example the recent studies on the reactivity of reactions confined inside micelles, vesicles and nanoreactors [START_REF] Walde | Enzymes inside lipid vesicles: preparation, reactivity and applications[END_REF][START_REF] Kim | Nanobiocatalysis and its potential applications[END_REF]). In ref. [START_REF] Chen | Diffusion into a nanoparticle with first-order surface reaction confined within a sphere[END_REF] the flux into the inner surface of the spherical cavity is seen to increase monotonically as the sink approaches the boundary of the cavity. The two limiting cases, concentric spheres and internally tangential, are associated with the minimum and maximum rate respectively. This result can be recovered easily by applying our method to calculate the reaction rate for a sphere Ω 1 of radius R 1 placed inside a spherical cavity Ω 0 of radius R 0 , as shown in Fig. 3.12. The boundary conditions read: If we let h 0 → ∞, this amounts to considering a fixed concentration of ligands at the membrane ∂Ω 0 .i.e. equal to the bulk concentration (see also the subsection 3.1.2). The system (3.49) in this case takes the form

∂u ∂ξ 1 -h 1 u ∂Ω 1 = 0 (3.56)
           -B l + ∞ n=l (l -h 1 ) (h 1 + l + 1) H n,l A n = 0 A l + l n=0 B n V n,l = δ l,0 (3.58)
with

V nl = (-1) l-n l n η l-n χ n+1 = (-1) l-n l n L R 0 l-n R 1 R 0 n+1 H nl = n l χ l η n-l (-1) n-l = n l R 1 R 0 l L R 0 n-l (-1) n-l ,
where, for the sake of clarity, we take the convention of placing the center of Ω 1 below the center of Ω 0 (see the scheme in Fig. 3.12). For the case of concentric spheres (L = 0 inside the system (3.58)), only one coordinate system is needed. This allows one to write the density function in Ω as:

u = A 0 + B 0 ξ 1 = A 0 + B 0 R 1 r .
Thus, in this case the monopole calculation coincides with the exact solution and the normalized rate can be explicitly computed by letting l = n = 0 in the system (3.58):

B 0 (1 + h 1 ) + h 1 A 0 = 0 A 0 + χB 0 = 1 (3.59) k k S = -B 0 = h 1 (1 + h 1 ) -h 1 χ . ( 3.60)
where k S = 4πDR 1 . For two totally absorbing boundaries (h 0 → ∞ and h 1 → ∞) the reaction rate is enhanced with respect to the unbounded domain (R 0 → ∞), and depends only on the ratio between the radii:

k k s = 1 1 -R 1 R 0 .
Before moving on to analyze more complex geometries it is instructive to examine another simple albeit non-trivial problem, namely the diffusive interaction between two sinks confined within a spherical cavity, separated by a permeable membrane from the outside (see figure 3.13). In the monopole approximation the system (3.49) with N = 2 and perfectly absorbing boundary conditions on each sphere reduces to:

     B (1) 0 + A 0 + B (2) 0 ǫ 2 = 0 B (2) 0 + A 0 + B (1) 0 ǫ 1 = 0 A 0 -B (1) 0 χ 1 -B (2) 0 χ 2 = 1, (3.61)
where χ i = R i /R 0 , ǫ i = R i /L. The normalized absorption rates to the two sinks Ω 1 and Ω 2 are:

           k 1 k S 1 = -B (1) 0 = a 1 1 1 + a 1 χ 2 + a 2 χ 1 k 2 k S 2 = -B (2) 0 = a 2 1 1 + a 1 χ 2 + a 2 χ 1 (3.62) with a 1 = 1 -ǫ 1 1 -ǫ 1 ǫ 2 , a 2 = 1 -ǫ 2 1 -ǫ 1 ǫ 2 .
The corresponding result in the unbounded domain (3.54) is recovered as R 0 → ∞.

To analyze the results, it is convenient to introduce a new function that measures the strength of the interaction between the two sinks: Here k [START_REF] Dhont | An introduction to dynamics of colloids[END_REF] (1+2) denotes the total rate into the two-sphere system, while k

H D = k [2] (1+2) k [1] S 1 + k [1] S 2 -1. (3.63)
[1]

S i = (1 -χ i ) -1
is the rate of a single sink inside the spherical cavity at the i-th position. All the rates are computed here in the monopole approximation. We will refer to the function H D as to the "diffusive interaction strength" because it quantifies the effects of the competition and cooperation between the reactive boundaries. The sign of H D identifies the nature of the interaction: a positive H D characterizes a cooperative interaction, the minus sign instead flags anticooperativity, while if H D = 0 there is no interaction at all. For two identical sinks of radius R, the previous expression gives:

H D = k [2] (1+2) 2k [1] S 1 -1 = a 1 1+2a 1 χ 1 1 1-χ 1 -1 = - R/L + 3R/R 0 1 + R/L + 2R/R 0 . (3.64)
highlighting that in the monopole approximation the behavior is always anticooperative. In computing the solution, the number of multipoles that we need to include increases as the ratio R/R 0 increases. As shown in figure (3.14), if the MOA is not a good approximation for this geometry, even for small sinks, the dipole approximation (obtained with N t = 2) already allows one to grasp the trend, at least for a reasonable range of the parameters. In figure (3.15) the ensemble of the two internal spheres is considered as a rigid body where the distance L between the centers is kept fixed and aligned to the z axis, while the two-body system is moved along the same line. The diffusive interaction is analyzed against as a function of distance between the origin and the center of the rigid body. The figure shows that the anticooperativity is maximum when the barycenter of the two-sink structure is placed in the origin, and it gets weaker as we move one of the sinks closer to the boundary (see again figure (3.13)).

Diffusive interactions and excluded-volume effects among multiple boundaries

The role of the environment's geometry on the kinetics of bimolecular reactions and the diffusive interaction between many reactive boundaries is a key aspect in modulating the rate constants of diffusion-influenced reactions. For this reason its effect has been widely studied in the literature, mainly experimentally and with numerical methods. The efforts for the minimization of the interaction of a biomolecule with secondary species, that can compromise the efficiency of the primary reaction, is a fundamental goal for example in drug development and industrial pharmacology. Thus many studies have been carried out to quantify the effects of geometry and reactivity characterizing the domain where the reaction takes place, mainly through numerical simulations of mathematical models. This investigation is the goal of a recent theoretical study [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF] that analyzes the screening of a sink, due to the presence of neighboring particles that can be either reactive or not with the same diffusing species. The analysis in Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF] is carried out through a numerical finite-element (FE) method, employed to solve the stationary diffusion equation in simplified geometries. This study provides clear-cut hints of the subtle effects brought about by the environment's geometry but also highlights the impossibility of brute-force numerical approaches to assess the impact of more crowded and sophisticated reactive environments. A spherical representation of the reactants alike to the one introduced in our study is used in Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF] to reproduce globular proteins and the reaction kinetics is described through a diffusion based model. The extremely high computational cost of the numerical FE procedure though, is a limit for the complexity allowed in the environment's representation. More specifically, the authors were able to analyze the diffusive interactions between multiple boundaries only for configurations involving at most 4 bodies, all of the same size. In the following, we will first show how our analytical theory is able to recover the numerical results reported in Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. Then we will generalize this problem to include more complex configurations comprising many spherical reactive boundaries. We will show that a comprehensive theory to quantify the role of the environment's reactive geometry can be built by unfolding this many-body problem through our theoretical method, since its generality allows one to compute the rate constant for a given reaction occurring in arbitrary reactive landscapes, made of multiple spherical boundaries of given size and reactivity in three dimensions. Moreover, we will demonstrate how ready-to-use closed approximation formulas can be derived easily in most cases, providing a valuable tool to capture the effects of competition and cooperation between neighboring particles.

Let us consider a ligand with bulk concentration c 0 that is absorbed by a sink of radius σ with its center at the origin. The screening action of N particles of specific reactivities h α can be cast in the form of a stationary problem for the normalized density of diffusing ligand molecules u = c(r)/c 0 :

               ▽ 2 u = 0 σ α ∂u ∂r α -h α u ∂Ωα = 0 ∀α = 1...N u| ∂Ω 0 = 0 lim r→∞ u = 1.
(3.65)

Even if our method could be employed to examine more complex reactive geometries, realized by assembling a large number of spherical boundaries of arbitrary size and reactivities, for the sake of clarity in comparing our results with the numerical data in [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF], we shall specialize here to the case of a central sink Ω 0 of radius σ surrounded by N identical spheres of radius σ 1 = λσ arranged randomly at a fixed distance d from the center of Ω 0 . Moreover, all the surrounding particles are taken to have the same reactivity h α = h ∀α. The solution can be written here as an expansion in series of irregular harmonics

u = 1 + N α=0 u - α (r α ), u - α = ∞ ℓ=0 ℓ m=-ℓ B α mℓ r ℓ+1 α Y mℓ (r α ). (3.66)
The expansion (3.66) guarantees that lim r→∞ u(r) = 1. Following the strategy introduced in section 3.1, we can claim that the boundary value problem is equivalent to the system:

-B α gq + (q -h α ) (h α + q + 1) δ g0 δ q0 + ∞ ℓ=0 ℓ m=-ℓ N β=0,β =α B β mℓ W αβgq mℓ = 0 (3.67)
for α = 0, 1, . . . , N , q = 1, 2, . . . , ∞ and g = -q, ..., q, where

W αβgq mℓ = (-1) q+g (ℓ -m + q + g)! (ℓ -m)!(q + g)! × R q α R ℓ+1 β L ℓ+q+1 βα Y m-g,ℓ+q (L βα ) (3.68)
with L βα = X β -X α , where X α is the position vector of the α-th screening sphere.

We recall that the normalized rate constant to the central sink is simply:

k 0 k S 0 = -B 0 00 .
In Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF] the solution is worked out by using the default direct linear solver FEniCS to compute the solution of the boundary problem (3.65) for N = 1, 2, 3, 4. Only the two limiting cases h α = 0 ∀α and h α → ∞ ∀α are considered, two choices that amount to consider, respectively, the screening spheres as perfectly reflecting and perfectly absorbing.

In figure 3.16 we compare the numerical FE calculation with the exact solution of (3.65) for N = 2, for both non-reactive and absorbing boundaries. It appears clear that the screening effect is harder to capture via a FE scheme in the case of reflecting obstacles than in the presence of competitive sinks. This discrepancy is in line with the general trend of the computational error involved in the numerical solution, which becomes larger as the reactivity of the neighboring spheres decreases The rate is plotted as a function of the normalized distance to the central particle. The analytical solution of (3.67) is compared to the numerical results taken from Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. The different symbols correspond to different angles between the neighbors. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. One remarks immediately that for every choice of the geometrical and chemical parameters the reaction rate to the central particle is reduced with respect to the flux to an isolated sink, because the external objects shield part of the ligand flux from it. It is possible to distinguish two main effects of the diffusive interaction. First of all, one has the excluded volume effect, which arises because the surrounding objects represent an obstacle that the diffusing molecules have to overcome in order to bind to the target site. However, a much stronger effect is represented by the competition among the three sinks for the common resources: the neighboring sites may bind the diffusing ligands, making them unavailable to the target receptor. The correction to the Smoluchowski rate appears to be short-range and non-dependent on the mutual positions of the neighbors for the case of non-reactive boundaries, while it becomes long-range and influenced by the mutual orientation in the case of reactive surrounding sites. The authors in Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF] suggested that, in the case of non-reactive neighbors, the rate function should depend linearly on the number of neighboring sites. By fitting the rate function with

k k S ∼ 1 -α(d)N,
where α(d) is a fit parameter dependent solely on the distance to the central sink, the authors of Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF] obtained for fixed choices of d a good agreement for the inert spheres, and a large discrepancy for the reacting boundaries. This observation is in line with the strong dependency on the configuration for the reactive neighbors case. This problem lends itself perfectly to illustrate the degree of analytical insight afforded by our theory. For this specific problem, it is expedient to look for solutions in the form of perturbative series.

k k S = ∞ j=0 ε j C j (N, θ, φ), (3.69)
where ε = σ/d is the perturbation parameter and, in principle, the coefficients C j (N, θ, φ) depend on the configuration of the surrounding spheres and are given by the sum of the many-body interactions. The perturbative procedure consists in expanding the variables of the linear system (3.67) in powers of ε, after truncating the multipole expansion to the value of N M which allows one to calculate the estimate (3.69) up to power j (see subsection 3.3.1 for detailed calculations). In the case of reflecting obstacles, one gets

k k S = 1 - λ 3 N 2 ε 4 - 2λ 5 N 3 ε 6 + . . . (3.70)
The approximation does not depend on the screening configuration and it is linear in N , as suggested in Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. However, in the subsequent discussion we will show how this only holds up to sixth order in ε: it can be seen from the expansion that the configuration enters explicitly for successive powers of ε.

A similar procedure in the case of N screening sinks yields

k k s = 1-λN ε+ λN +λ 2 N α,β=1 β =α 1 Γ αβ ε 2 -λ 2 N 2 +λ 2 N α,β=1 β =α 1 Γ αβ +λ 3 N α,β,δ=1 β,δ =α 1 Γ αβ Γ αδ ε 3 +. . . (3.71) where Γ αβ = 2 sin(ω αβ /2),
ω αβ being the angle formed by the sinks α and β with respect to the origin. Eq. (3.71) makes it very clear that the configuration of competitive reactive boundaries does influence the screening effect on the central sink. In fact, in this case the screening boundaries compete with each other for the ligand flux, and the mutual orientation of the neighbors has to be taken into consideration even for small values of ε, because the configuration enters the picture already at the second order term of the expansion through the functions Γ αβ . The second-order term in the expansion shows clearly the competition effects between the neighboring reactive particles: the closer are the bodies (the smaller is Γ(α, β)), the stronger is the positive second-order contribution to the correction. Thus, the more the neighbors are packed closely, the more the absorbing rate of the central sink is increased. This is evident in Fig (3.16): the screening effect on the central sink is stronger when the angle between the neighbors is larger. In this sense the competition among neighboring reactive particles can prove beneficial to an observed site. A clear Exact solution Linear 2nd order 3rd order signature of this competitive effect is also that the corrections in Eq. (3.71) alternate in sign. This observation sheds considerable light onto the many-body character of the rate constant, whose perturbative series is alternatively reduced by the diffusive interactions between the screening boundaries and the sink (shielding ligand flux from it) and increased by the diffusive interactions among the screening particles (shielding flux from each other, see Fig. 3.17). On the contrary, the screening action of inert obstacles is largely dominated by the excluded-volume effect, and thus can only yield negative corrections at all orders.

Due to its perturbative nature, eq. (3.71) can be used to quantify the shielding action of specific 3D arrangements of sinks only for N ε ∝ N/d ≪ 1. However, it still provides a powerful analytical tool to compare different geometries, as the perturbative rate is always proportional to the true rate. In Fig. 3.18 we compare the perturbative expansion (3.71) of the rate constant corresponding to given configurations with the exact value. More specifically we plot the rate constant for a sink surrounded by competitive perfectly absorbing spheres, computed through eq. (3.71) truncated at the second and third order, as a function of the exact results, obtained by solving the linear system. The purpose of these graphs is to show that even if the perturbative approximation generally overestimates the exact results, unless N ε ∝ N/d ≪ 1, it yields a result which is to an excellent extent proportional to it. Thus, eq. ( 3.71) provides an effective figure of merit to compare the effect of different conformations, e.g. with the aim of designing the special configurations that minimize or maximize the screening effect on the central sink for given values of N and d.

The competitive and excluded volume effects are highlighted in Fig. 3. [START_REF] Fanelli | Diffusion in a crowded environment[END_REF], where the influence of the configurations is shown and the calculation is compared to the perturbative approximation, for the case of reflecting bodies, and to the MOA approximation of the system (3.67) averaged over many different configurations of the reactive neighbors. The average shielding action exerted by N equidistant reactive sinks can be in fact captured analytically in the monopole approximation, i.e. by keeping only the ℓ = 0 and q = 0 terms in eqs. (3.67). The ensuing reduced system can be averaged over different configurations in the hypothesis of vanishing manybody spatial correlations, i.e. by averaging the equations over the multi-variate totally uncorrelated probability density of the angles between each pair of spherical sinks The result for the mean rate is (see subsection 3.3.2 for the details) 3.73) provides an very good estimate of the configurational averages of the exact results at separations greater than a few diameters, highlighting the dramatic screening action of competitive reactive boundaries with respect to inert obstacles. Furthermore, a simple analysis of the rate fluctuations over the configuration ensembles at fixed d allows one to gauge how sensitive the competitive screening is to the specific 3D arrangement of the sinks. Remarkably, this analysis reveals stretches between the minimum and maximum rates for a given value of d as high as 50 % of the average (see shaded bands in Fig. 2). More precisely, we remark that the variability associated with different geometries is greater (i) at short distances and (ii) for few screening particles.

P( - → w ) ≡ α =β P (ω αβ ), ( 3
k k S = 1 -λε[N -(N -1)(1 -λε)] 1 -λε[N ε -(N -1)(1 -λε)] . ( 3
To stress the importance of the configuration and to relate it to the number of the particles, it is interesting to extend the previous analysis to the case of N spheres randomly arranged within a spherical volume around the central sink in a shell centered on the selected sink, and to compare the absorbing rates obtained by varying the packing fraction and the configuration of the neighbors. Fig. 3.20 shows in a clear fashion that the packing fraction alone cannot explain the modulation of the reaction rate: the configuration of the spheres in the volume does matter. Interestingly, allowing the particles to occupy an arbitrary position in a volume around the sink, this effect is evident for both inert and reactive objects. The large variance in the results obtained by varying the arrangement of the agents appears related to the average distance of the crowders from the center as well as to the homogeneity in their distribution. By choosing the standard deviation between the N surrounding spheres as a measure of order for the crowders' assemblies, we see that the heterogeneity in the configuration correlates positively with the modulation of the rate, once the packing fraction is fixed. The normalized rate into the target sink is found to increase for less ordered configurations of the screening spheres, both in the case of reflecting and absorbing neighbors. We can explain this correlation with two effects: as the arrangements become more regular, the increased probability of occurrence of holes in the screening configurations facilitates the diffusion of the ligands. In addition to this effect, the higher competition between the spherical boundaries increases the efficiency of the target sink. The latter effect appears to be stronger, since the modulation of the rate is more evident when we consider the case of screening sinks (right panel in Fig. 3.20). Even if the neighbors shield almost all the flux from the sink (we are considering here N = 34 reacting boundaries), the relative difference between the rates obtained in the less and most ordered settings is as high as 80%.

Outline and details of the perturbative expansions

In order to explain the procedure followed to derive eqs. (3.70) and (3.71), we will describe here in some detail the steps that led us to compute the first two nontrivial terms for the non-reactive boundaries expansion. The extension to higher orders involves longer but straightforward calculations along similar lines. Let us rewrite the linear system obtained for the case of N reflecting spheres, and let us single out the coefficients of the central sink B 0 gq from the B α gq that characterize the reflecting boundaries, following expressions for all the coefficients of the linear system:

             -B 0 gq -δ (q,g)=0,0 + ∞ n=1 n m=-n N β=1 B β mn W 0βgq m,n = 0 -B α gq + q q + 1 ∞ n=0 n m=-n B 0 mn W α0gq m,n + ∞ n=1 n m=-n N β=1,β =α B β mn W αβgq m,n = 0 (3.74) where W αβgq m,n = (-1) q+g (n -m + q + g)! (n -m)!(q + g)! σ q α σ n+1 β L n+q+1 βα Y m-g,n+q (L βα ) (3.
B 0 m,n = ∞ j=0 ε j Q 0,j m,n B α m,n = ∞ j=0 ε j Q α,j m,n (3.76)
where Q α,j m,n is the j-th coefficient of the expansion. We have to take into account the number of multipoles which is necessary to obtain the desired order of approximation. The fourth order of the expansion in powers of ε requires only the terms

= λ 3 N 2 (3.78)
The corresponding expression for a collection of N reactive boundaries can be obtained with the same procedure, the only additional difficulty being the contribution of the pairwise distance between the neighboring boundaries L αβ , which enter the expansion already in the first-order terms. To carry out the approximation, it is necessary to express L αβ as a function of the perturbative parameter by using basic trigonometry, i.e.

σ L αβ = ε 2 sin ω αβ 2 (3.79)
where ω αβ is the angle formed by the sinks α and β with respect to the central sink.

As a consequence, the perturbative approximation will depend on the configuration of the surrounding boundaries.

Our perturbative procedure has shown that, in case of N reflecting boundaries, the rate constant is linear in N and depends only on the distance between the obstacles and the sink. The configuration does not appear explicitly. However, this is only the case at large enough separations, as the specific three-dimensional arrangement of obstacles indeed enters the perturbative expansion at the 7-th order in ε = σ/d. In order to prove this statement, let us consider two distinct planar configurations of the reflecting spheres, lying at the vertices of an equilateral triangle (configuration C 3 ) and of a square (configuration C 4 ). The seventh-order correction for a collection of spheres lying in the z = 0 plane, θ = π/2, reads

Q 0,7 00 = - 1 16 N β=1 N β=1,β =γ 1 Γ(β, γ) 3 e i(φ 0β -φ 0γ ) + e -i(φ 0β -φ 0γ ) +3 e i(φ 0β +φ 0γ -2φ γβ ) + e -i(φ 0β +φ 0γ -2φ γβ ) (3.80) = - 1 8 N β=1 N β=1,β =γ 1 Γ(β, γ) 3 cos(φ 0β -φ 0γ ) + 3 cos(φ 0β + φ 0γ -2φ γβ ) .
One immediately recognizes from eq. (3.80) that the configuration enters the picture through the parameter Γ(β, γ) and through the azimuthal angles φ αβ . Evaluating the corrections for the two regular polygons, Q 0,7 00 (C 3 ) and Q 0,7 00 (C 4 ), we find

Q 0,7 00 (C 3 ) N = 7 24 √ 3 = Q 0,7 00 (C 4 ) N = 3 2 √ 2 + 1 4 .
This proves that for orders higher than six in ε, the perturbative expansion of the rate constant to the sink is no longer independent of the specific arrangement of reflecting obstacles. Nevertheless, for a wide range of the parameters ε and N , one can neglect the dependency on the configuration for the case of reflecting boundaries. We have shown that this is no longer true for a collection of reactive neighbors. In fact, already the coefficient of the second order term in the expansion (3.71) do not have a linear dependence on the number of neighboring boundaries, unlike the corresponding formula for the reflecting spheres. To justify this claim, we can estimate explicitly the leading power in N of the second-order coefficient for a particular domain. For the sake of simplicity, we will confine ourselves to a planar configuration with the N centers of the spheres lying at the vertices of a regular polygon and also take λ = 1. In this case, the pairwise distances for each pair Ω α , Ω β can be written explicitly, so that

k k s = 1 -N ǫ + N + N α=1 N β =α=1 1 Γ(α, β) ǫ 2 + o ǫ 2 = 1 -N ǫ + N + N N -1 k=1 1 2 sin(πk/N ) a(N ) ǫ 2 + o ǫ 2 . (3.81)
It can immediately be concluded that the order in N of Q 0,2 00 = N + a(N ) is higher than two. One has lim

N →∞ a(N ) N 2 = lim N →∞ 1 N N -1 k=1 1 2 sin(πk/N ) = 1 0 dx 2 sin(πx) = π 2 0 dx sin(x) ≥ π 2 0 dx x → ∞
A more precise estimate can be obtained using the approximation

1 sin(πk/N ) = N πk + r(k, N )
where r(k, N ) is the rest. Assuming N odd and using the properties of the harmonic series (for large values of N one has N k=1 k -1 ∼ ln N ), we get

N N -1 k=1 1 2 sin(πk/N ) = N N -1 2 k=1 1 sin(πk/N ) ∼ N 2 π (N -1)/2 k=1 1 k ∼ N 2 π ln N -1 2 .
Since it is possible to show that N (N -1)/2 k=1 r(k, N ) has order N 2 , we conclude that, at least in the planar ordered configurations considered here, the second coefficient of the expansion has order N 2 ln N . In fact, we have

(N -1)/2 k=1 r(k, N ) = (N -1)/2 k=1 1 sin(πk/N ) - N πk and lim N →∞ 1 2N N k=1 1 sin(πk/2N ) - 2N πk = 1 π π 2 0 1 sin(y) - 1 y dy < ∞ because lim y→0
ysin y y sin y = 0.

Monopole approximation

In this subsection we will show how to compute the monopole approximation. This is in general a very powerful approximation, especially when the reactive landscape is dominated by a large number of sinks. Here we will illustrate this technique with reference to the derivation of eq. (3.73). Let us write explicitly the linear system ∀α = 0, . . . N, ∀q = 0, . . . ∞, ∀g = -q . . . q:

-B α gq + (q -h α ) (h α + q + 1) δ g0 δ q0 + ∞ ℓ=0 ℓ m=-ℓ N β=0,β =α B β mℓ W αβgq mℓ = 0
for a configuration of N sinks randomly placed at a distance d from the central sink in the monopole approximation, that is, ℓ = q = 0. We get ∀α = 1, 2, . . . , N

             B 0 00 = -1 - N i=1 B i 00 W 0i00 00 B α 00 = -1 -B 00 W α000 00 - N β =i=1 B β 00 W αβ00 00 ∀α = 1 . . . N.
We recall that in general the rate to the α-th sink is:

k α k Sα = -B α 00 .
After summing the second equation over i, and letting x = B 0 00 and y

= N i=1 B i 00 we obtain        x + 1 + λǫy = 0 y + N + N ǫx + λǫ N i=1 N β =i=1 B β 00 Γ α,β = 0.
Physically, the variable x is the negative of the rate constant to the central sink, k/k S = -x, while the variable y stands for the total rate constant for the ensemble of N screening spheres, namely

N α=1 k α k Sα = -y.
Obviously the solution will depend on the configuration of the spherical boundaries, which is embodied in the functions Γ α,β . The value of the mean rate can be obtained by averaging both equations over the multi-variate probability density (3.72) of the angles ω αβ between each pair of spherical sinks. The excluded-volume constraint among the screening sinks requires two given particles to lie at a distance greater than or equal to their diameter. This means that 2 arcsin(λε) ≤ ω αβ ≤ π ∀α, β, with λ = σ 1 /σ and ε = σ/d. The average of the configuration-dependent terms gives

1 Γ αβ = π 2 arcsin λε P (ω) dω 2 sin(ω/2) = 1 -λε.
Noting that

N i=1 N β =i=1 B β 00 = (N -1) N i=1 B i 00 = (N -1)y we finally get k k S = -x = 1 -λε[N -(N -1)(1 -λε)] 1 -λε[N ε -(N -1)(1 -λε)] . (3.82)
This expression reduces to a well-known formula for a system of two identical sinks [START_REF] Deutch | Competitive effects in diffusioncontrolled reactions[END_REF], k/k S = 1/(1 + ε).

In figure 3.19 the averaged MOA approximation is compared to the averaged exact results for reactive sinks in the particular case of λ = 1, which is the setting analyzed in Ref. [START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. The picture clearly shows how, for the choice σ = σ 1 , the rate to the central particle increases monotonically with its distance from the screening spheres.

It is instructive to analyze to which extent this behavior changes if we allow the central sink to have a different size from the screening boundaries. It is straightforward to prove that the averaged MOA displays a minimum as a function of the distance d for certain choices of the parameters N and λ. Since ε ∝ d -1 , the condition for the existence of a stationary point as the distance d is changed reads

d dε k k S = - λN (ǫ(λǫ(λ(N -1) + 1) -2) + 1) (λǫ(λ(N -1)ǫ + N (ǫ -1) + 1) -1) 2 = 0. (3.83)
From eq. (3.83) it follows that a minimum exists for

λ 2 (N -1) + λ -1 ≤ 0. (3.84)
In this case, the distance d * at which a minimum is found is

d * = σ 1 -λ[1 + λ(N -1)] + 1 . (3.85)
If a minimum exists, it has to occur at a distance greater than or equal to the contact distance between the central sink and the screening particles. Thus, one should enforce the condition d * ≥ σ + σ 1 . Therefore, in view of expression (3.85), it follows that one has to complement eq. (3.84) with the additional requirement

1 -λ[1 + λ(N -1)] + 1 ≥ 1 + λ. (3.86)
The solution to the system of inequalities (3.84) and (3.86) reads

λ ≤ λ * (N ) ≡ ( √ 4N + 1 -1)/(2N ) < 1 for fixed N N ≤ N * (λ) ≡ (1 -λ)/λ 2 for fixed λ (3.87)
The non monotonicity is another non-trivial effect of the mutual screening and of the diffusive interaction among the reactive boundaries (see the discussion in section 3.2.1). This effect is shown in Fig. 3.22, where the exact rate to a sink surrounded by smaller reactive neighbors is compared with the MOA approximation for an average configuration of the particle. We can note that the non monotonicity is here more evident than in the two-sink configuration considered in 3.2.1. Upon increasing the distance d, in fact, the total flux to the many screening sinks will increase because their active surfaces get larger and they also get farther apart from each other. In the following chapters 4 and 5 we will present the application of our model to two examples of diffusion influenced reactions. In chapter 4 we aim at characterizing, as a function of the structural parameters, the efficiency of complex architectures of core-shell nanoreactors where many catalytic metallic nanoparticles are embedded in a polymeric material around a plastic core to form a complex reactor. Chapter 5 deals with the biological, naturally occurring binding reaction between an antigen and an antibody, and more generally points to a way to adapt our procedure to account for ligands binding to complex non-spherical large biomolecules. The approach for modeling the antigen-antibody interaction combines a coarse grained dynamical representation of the protein with the definition of the specific multiple boundaries problem.

Chapter 4

Catalytic rate constant of complex core-shell nanoreactors A thriving field in modern nano-sciences is the design and fabrication of novel composite catalytic structures. In particular, several strategies have been devised to employ large number of catalytic metal nanoparticles to assist a wide array of oxidation, reduction and hydrogenation reactions. The interest in composite organicmetal active nanostructures arises due to the fact that metals often have totally different properties at the nano-scale [START_REF] Herves | Catalysis by metallic nanoparticles in aqueous solution: model reactions[END_REF]. Several metals have been studied and adapted to be used as nanocatalysts: gold, silver, palladium, silica particles have been assembled together within organic meshes of different sorts in various configurations and many efforts have been made to optimize the structure of these composite catalytic systems. The catalytic efficiency of an assembly of metallic nano-particles is directly connected to their global surface area, because a larger global surface allows more reactions to occur at the same time. However, the resulting need for assembling a large number of nanoparticles poses a problem, because the nanoparticles in liquid solutions tend to aggregate at large concentrations. To keep them separated in a fixed position, the catalysts can be either endowed with protective layers or fixed on the surface of static objects, or even immobilized in colloidal carrier systems. A recent innovative line of technology is to use organic meshes whose properties are tunable in response to the variation of external parameters, such as PH and temperature. In this chapter we will focus on particular reactors consisting of several gold, silver or palladium nanoparticles immobilized in a microgel made of poly-N-isopropylamcrylamide (PNIPAM), a polyelectrolite which can shrink and swell following changes in temperature. The microgel is fixed around a spherical plastic core, thus forming a shell where the catalysts are embedded and where the reactants diffuse until they get in contact with the nanoparticles (see the scheme in Fig. 4.1). In addition to keep the particles fixed throughout the reaction, the use of a termosensitive colloid provides a powerful way of controlling the reaction rates, since the varying thickness and the modified chemical properties of the carrier influence the mobility of substrate molecule within the composite structure.. In particular, the PNIPAM microgel undergoes a phase transition at 32 • C: increasing the temperature past this critical value, the gel passes from a swallen to a shrunk state, by expelling water. This property has a complex effect on the reaction kinetics: on one the hand the diffusion coefficient is decreased, because the medium becomes thicker. On the other hand, depending on the specific reactant and its affinity with water, the change of phase can have a positive or negative effect on the reactant mobility. More precisely, the solvation free energy of substrate molecules has been shown to play a crucial role in controlling the overall catalytic activity of the nano-reactor [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]. In reason to the ever growing number of applications of these technologies in chemical and biological applications, the need for a deep understanding of how the structural parameters of these composite catalytic scaffolds influence the efficiency of the catalysis has recently become a topic of great interest. The analythical tool that we have developed in our work provides the perfect instrument to characterize the best setting in terms of the dimension and arrangement of the components, temperature range and choices of the materials. More specifically, we set out to employ the model described in the previous sections to compute the overall reaction rate constant of a composite core-shell nanoreactor, consisting of a PS core surrounded by a hydrogel layer of varying thickness, where a given number of small reactive spheres are embedded at prescribed positions (see the sketch in figure 4.1). These are endowed with a prescribed intrinsic rate constant. Both the diffusion coefficient in the hydrogel and in the bulk enter our model explicitly. As a result of our analysis we not only provide a way of computing the reaction rates corresponding to various configurations exactly. We also worked out a simple but very accurate analytical formula for the overall rate constant, derived within the monopole approximation. Our analytical treatment provides a very useful and flexible tool for the design of optimally performing composite nanoreactors with different polymer shells and carrying different kinds of nano-catalysts.

Mathematical model for the catalytic reaction

We model a nanoreactor consisting of a polysterene (PS) core surrounded by a microgel (or polymer brush) layer as two concentric spheres centered at the origin of a 3D frame, as depicted in Fig. 4.2. We denote with R S and R 0 the core and shell radius, respectively. The shell carries N small nano-catalysts (metal nanoparticles or enzymes) that we model as spheres of radius a. For the sake of simplicity, we label the PS core as the inner sphere with position vector L 1 = 0 while the vectors L α , α = 2, 3, . . . , N + 1 denote the position of the N nanocatalysts. We want to compute the reaction rate constant for reactions taking place at the surface of the catalyst spheres. These are endowed with an intrinsic rate constant k * , which in principle could be a function of temperature of the Arrhenius-type [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]:

k * ∝ exp - ∆E k B T
where ∆E is the activation energy required to trigger the reaction. Let us denote with S 0 ≡ {r 0 , θ 0 , ϕ 0 } the reference frame with the origin at the nanoreactor center and with S α ≡ {r α , θ α , ϕ α } the N reference frames with the origins at the nanospheres centers and the axes parallel to S 0 . This formally defines the following 3D domains

Ω + = {r 0 ∈ (0, R 0 ], θ 0 ∈ (0, π), ϕ 0 ∈ (0, 2π)} \ ∪ α Ω α Ω -= {r 0 ∈ [R 0 , ∞), θ 0 ∈ (0, π), ϕ 0 ∈ (0, 2π)} (4.1)
where Let ρ B denote the bulk density of reactants and let us introduce the timedependent normalized density u(r, t) = ρ(r, t)/ρ B . As usual, we assume that the system relaxation time for the diffusive flux of B particles (the ligands) t D ≃ (R 0 -R S ) 2 /D i is small enough to neglect time-dependent effects. Hence, in the absence of external forces, the diffusion of ligands with normalized number density u(r) is described by the steady-state diffusion equation

Ω 1 = {|r 0 | < R S }
∇ • [D(r)∇u(r)] = 0 in Ω = Ω + ∪ Ω - (4.2)
with

D(r) = D i in Ω + (microgel) D o in Ω -(bulk) (4.3)
and with the customary bulk boundary condition

lim |r|→∞ u(r) = 1. (4.4)
The difference with the setting considered up to now is that here the problem is defined in both the internal (confined) Ω + and external (unbounded) Ω -region.

It is well known from the general theory of partial differential equations that the classical solution (twice continuously differentiable in Ω and continuous on Ω) of the stationary diffusion equation (4.2) does not exist in the whole domain Ω [START_REF] Ladyzhenskaya | Linear and quasilinear elliptic equations[END_REF]. Therefore one should consider the function

u(r) = u + (r) for Ω + (microgel) u -(r) for Ω -(bulk). (4.5)
Accordingly, we should impose a condition for the substrate concentration field at the bulk/microgel interface, ∂Ω 0 ≡ {r 0 = R 0 }. It has been demonstrated recently that a key factor controlling the overall reaction rate is the solvation free-energy change ∆G, a quantity that describes the partitioning of the reactant in the microgel versus bulk [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]. It is possible to show in the case of a single catalyst that a free-energy jump at the solvent-microgel interface can be accounted for by a modified ligand density in the microgel versus bulk, ρ B → ρ B exp(-β∆G) [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]. Here we assume that such description holds unchanged in the case of multiple catalysts. Accordingly, we require u + -λu - ∂Ω 0 = 0 (

where λ = exp(-β∆G). (

Furthermore, the following continuity condition for the local diffusion fluxes should also hold at the bulk/microgel interface

∂u - ∂r 0 -ζ ∂u + ∂r 0 ∂Ω 0 = 0 (4.8)
where we have introduced the diffusion anisotropy parameter

ζ = D i D o (4.9)
Conditions (4.6) and (4.8) for D i = D o are often called the weak discontinuity conditions for the concentration field u(r). Finally, reflecting boundary conditions should hold at the surface of the PS core, i.e.

∂u + ∂r 0 r 0 =R S = 0. (

We are interested in describing the pseudo-first-order irreversible diffusion-influenced reaction between the N nano-catalysts C encapsulated in the microgel and the reactants B which are freely diffusing in the bulk and in the microgel. The process can be cast in the usual scheme:

C + B k D ---⇀ ↽ --- k -D CB * k * -→ C + P (4.11)
where as usual CB * denotes the encounter complex, k D and k -D are the association and dissociation diffusive rate constants, respectively, and k * is the intrinsic rate constant of the chemical reaction occurring at the sufrace of the metal nanoparticles.

In substrate excess conditions, the reaction (4.11) is equivalent to the following pseudo-first-order kinetic scheme

C + B k -→ C + P. (4.12)
Our final goal is to compute the global diffusion-influenced rate constant k which we can calculate as

k = N +1 α=2 ∂Ωα D i ∂u + ∂r α ∂Ωα dS (4.13)
and to characterize its dependence on the physical-chemical parameters of the nanoreactor. We stress that the introduced schematization of the problem holds under the excess reactant condition c ≪ ρ B , where we denote with c the concentration of the catalysts. Introducing in the following the usual dimensionless spatial variables ξ 0 = r 0 /R 0 , ξ 1 = r 0 /R S and ξ α = r α /a for α = 2, 3, . . . , N + 1, our problem can be cast in the following form

∇ 2 u ± = 0 in Ω ± (4.14a) ∂u + ∂ξ α -hu + ∂Ωα = 0 α = 2, 3, . . . , N + 1 (4.14b) lim ξ 0 →∞ u -(ξ 0 ) = 1 (4.14c) ∂u + ∂ξ 1 ξ 1 =1 = 0 (4.14d) u + -λu - ∂Ω 0 = 0 (4.14e) ζ ∂u + ∂ξ 0 - ∂u - ∂ξ 0 ∂Ω 0 = 0 (4.14f)
The parameter

h = k * 4πaD i ≡ k * k + S , (4.15)
as in the previous sections is meant to gauge the character of the reaction, and it is normalized to the Smoluchowski rate constant of a single nanosphere embedded in the microgel, k + S = 4πaD i . We look for the solutions of the above problem in the bulk and in the microgel as linear combinations of regular and irregular harmonics:

u + (r) = ∞ n=0 n m=-n A mn ξ n 0 Y mn (r 0 ) + N +1 α=1 ∞ n=0 n m=-n B α mn ξ -n-1 α Y mn (r α ) (4.16a) u -(r) = 1 + ∞ n=0 n m=-n E mn ξ -n-1 0 Y mn (r 0 ) (4.16b)
where ξ 0 = r 0 /R 0 , ξ 1 = r 0 /R S , ξ α = r α /a for α = 2, 3, . . . , N and A mn , B α mn and E mn are N + 3 infinite-dimensional sets of unknown coefficients.

In order to determine the unknown coefficients we express the solution in the local coordinates on every boundary (the N + 1 spherical surfaces ∂Ω α ) and at the microgel-bulk interface ∂Ω 0 , where we impose the weak discontinuity conditions for the ligand density field.

With the same strategy introduced in section 3.1, and eliminating the coefficient E mn , we obtain the following linear equations:

1 λ + q q + 1 A gq + 1 λ -ζ N β=1 q n=0 n m=-n B β mn V β,m,n g,q I {g-(q-n)≤m≤g+(q-n)} = δ g0 δ q0 (4.17a)
-B α gq + (q -h α ) (h α + q + 1) ∞ n=0 n m=-n A mn H (α,g,q) m,n I q≤n + N β=1,β =α B β mn W (α,β,g,q) m,n = 0, (4.17b)
where h 1 = 0, h α = h for α > 1. Eqs. (4.17) hold ∀ q ∈ [0, ∞) with α = 1, 2, . . . , N +1 and g = -q, -q + 1, . . . , q -1, q. The matrices V, H, W read

V α,m,n g,q = (-1) q-n+m-g (q -g)! (n -m)!(q -n + m -g)! η q-n 0α χ n+1 α Y m-g,q-n (-L α ) (4.18a) H (α,g,q) m,n = n + m q + g χ q α η n-q 0α Y m-g,n-q (L α ) (4.18b) W (α,β,g,q) m,n = (-1) q+g (n -m + q + g)! (n -m)!(q + g)! η -(n+q)-1 βα χ q α χ n+1 β Y m-g,n+q (L βα ) (4.18c)
where L αβ = L β -L α (according to this notation L 0α = L α ) and

η αβ = η βα = L αβ R 0 χ α = R α R 0 . (4.19)
Here, for the sake of trasparency, we pose R 1 = R S (radius of the PS core) and R α = a, for α > 1 (radius of the nanocatalysts).

If the multipole expansions are truncated at N M multipoles, the system (4.17)) comprises (N +2)(N M +1) 2 equations, which can be easily solved numerically. Once the coefficients have been determined, the overall rate constant can be obtained easily as:

k = N α=2 k a = -k + S N +1 α=2 B α 00 . (4.20)
The system one has to solve has the following structure

              1 λ + ζq q + 1 1 1 λ -ζ V 1 1 λ -ζ V 2 . . . 1 λ -ζ V N +1 H 1 -1 W 1,2 . . . W 1,N +1 H 2 W 2,1 -1 . . . W 2,N +1 . . . . . . . . . . . . . . . H N +1 W N +1,1 W N +1,2 . . . -1               ×                       A 00 . . . A N M N M B 1 00 . . . B 1 N M N M . . . B N +1 00 . . . B N +1 N M N M                       =                       1 . . . 0 0 . . . 0 . . . 0 . . . 0                      
The above described procedure allows one to compute the reaction rate constant for a given geometry of the system and for given values of the relevant physico-chemical parameters, with arbitrary precision. The monopole approximation of the solution of the system for a given configuration can be obtained by truncating the expansion to q = n = 0. The ensuing equations read

1 λ A 00 + 1 λ -ζ N β=1 B β 00 V β00 00 = 1 (4.21a) -B α 00 - (h α ) (h α + 1)
A 00 H (α00) 00

+ N β=1,β =α B β 00 W (αβ00) 00 = 0 α = 1, 2, . . . , N (4.21b)
Recalling the definitions of the matrices, we have

V β00 00 = a/R 0 , H (β00) 00
= 1, W (αβ00) 00 = a/L βα , and ζ = D i /D o , so that the system takes the form

                   1 λ A 00 + 1 λ -ζ a R 0 N β=1 B β 00 = 1 B α 00 + h α 1 + h α A 00 + N β=1,β =α B β 00 a L αβ = 0 α = 1, 2, . . . , N . (4.22)
Since B α 00 = -k α /k + S , the overall rate constant of the nanoreactor can be computed simply as k = -N +1 β=1 k S B β 00 (note that B 1 00 = 0 is identically zero as the PS core is modeled as a reflecting sphere). Moreover, we can average the system (4.22) over the nanoparticles configurations P(L 2 , L 3 , . . . , L N +1 ). We therefore get from Eqs. (4.22)

         1 λ A 00 - 1 λ -ζ a R 0 k k + S = 1 k k + S - h 1 + h N A 00 -(N -1) k k + S a L αβ = 0 (4.23)
where we have taken h α = h = k * /k + S ∀ α as the N catalysts are identical. We note that a more useful quantity is the total rate constant normalized to the Smoluchowski rate of an isolated sink of the same size as the whole nanoreactor in the bulk, k - S = 4πD o R 0 . Eqs. (4.23) can be solved straightforwardly, which gives

k k - S = N k * a R 0 ζe -β∆G k + S + k * 1 + (N -1) a L αβ - N a R 0 1 -ζe -β∆G (4.24)
where we recall that λ = e -β∆G , ζ = D i /D o and k + S = 4πD i a. The quantity a/L αβ represents the average inverse inter-catalyst separation, which can be computed analytically under the reasonable assumption that spatial correlations in the catalysts configurations are negligible. The probability density for the position of the centers of any two particles in a shell of radii R S + a and R 0a is thus given by

P(r, ρ) = 3 (R 0 -a) 3 -(R S + a) 3 2 r 2 ρ 2 sin θ 4π (4.25)
so that the average inverse inter-catalyst separation is

a L αβ = 9 a 2[(R 0 -a) 3 -(R S + a) 3 ] 2 (4.26) × R 0 -a R S +a r 2 dr R 0 -a R S +a ρ 2 dρ π 0 sin θ r 2 + ρ 2 -2rρ cos θ dθ = 2(1 -ε) 5 -5(1 -ε) 2 (γ + ε) 3 + 3(γ + ε) 5 (1 -ε) 6 -2(1 -ε) 3 (γ + ε) 3 + (γ + ε) 6 3a 5R 0 : = ε C(ε, γ) (4.27)
where γ = R S /R 0 denotes the fraction of the nanoreactor size occupied by the PS core and ε = a/R 0 is the non-dimensional size of each catalyst. We see that, since ε ≪ 1, one has 1 + ε/3 C 6(1 + ε)/5, i.e., C is of the order of unity, 1.005 C 1.217 (taking ε ≈ 0.0146 from experiments [START_REF] Mei | Catalytic activity of palladium nanoparticles encapsulated in spherical polyelectrolyte brushes and core-shell microgels[END_REF]). This has the interesting implication that the exact size of the core does not play a big role for the overall rate process for relevant values of the physical parameters, as we discuss further below.

If the catalytic action exerted by the metal nanoparticles encapsulated in the microgel is fast with respect to diffusion, i.e. k * ≫ k + S , expression (4.24) can be simplified by taking the limit k * → ∞. This yields the expression for the fully diffusion-controlled rate

k k - S = N ε ζe -β∆G 1 + (N -1)ε C(ε, γ) -N ε 1 -ζe -β∆G . ( 4.28)
Eq. (4.24) is a very important result, as it allows one to compute the overall reaction rate as a function of all the relevant structural and physico-chemical parameters that control the performance of the nanoreactor. We stress that in this formula the ligand diffusion coefficient in the microgel, D i , the radius of the nanoreactor, R 0 , and the effective intrinsic rate constant, k * , are all functions of temperature. In the next section we will discuss more in depth the properties and the usefulness of the analytical solution.

Discussion of the analytical solution

We now discuss the essential features of the diffusion-controlled rate in Eq. (4.28). For a random distribution of the nanoparticles in the shell we found a/L αβ = Ca/R 0 , where C depends on the ratio between the radius of the PS core and that of the nanoreactor and is of the order of unity. Hence, we find that the exact size of the core actually does not play a significant role for the diffusion-controlled rate for relevant values of the physical parameters, i.e., (weak) attraction to the hydrogel ∆G < 0 and decreased internal diffusion ζ < 1. The most pronounced effect of the core corresponds to the case ζ = 1 (bulk diffusion in the shell) and ∆G = 0. Then, comparing the limits of vanishing core, (γ = 0 which yields C(γ = 0) ≃ 6/5) to the thin-shell limit, (γ ≃ 1 which yields C(γ = 1) ≃ 1), shows that the two corresponding values of the overall maximum rate k/k - S (for N → ∞) vary from 5/6 to 1, i.e. they differ by only 20 %. In other words, if the nanoparticles are distributed uniformly inside the nanoreactor or only in the surface shell region makes only a difference of at most 20% in the rate. In the following we neglect these effects, as we move to a qualitative discussion of the more significant effects.

As we see from Eq. (4.29) the maximum achievable rate is k = k - S = 4πD 0 R 0 , that is, the Smoluchowski rate of a sink of size equal to that of the total nanoreactor, i.e. the nanoreactor should be big for high activity. In the limit of small nanoparticles-to-nanoreactor size ratio, ε ≪ 1, Eq. (4.28) can be simplified to the following form

k k - S = N ε ζe -β∆G 1 + N ε ζe -β∆G . (4.29)
Let us recall the important parameters, that is, the nanoparticles-to-nanoreactor size ratio 0 < ε = a/R 0 ≪ 1, the number of nanoparticles N , the scaled reactant mobility inside the shell 0 < ζ = D i /D 0 1, and finally the change in solvation free energy ∆G for the reactants upon entering the hydrogel. Clearly, if the mobility vanishes, ζ ≪ 1 or the solvation free energy change ∆G ≫ k B T is substantially repulsive, the reaction is significantly slowed down. However, in realistic systems the mobility will be certainly slowed down to some extent but not vanish. ∆G may be even negative (attractive) if the reactant interacts favorably with the polymer as found for rather hydrophobic reactants and collapsed PNIPAM-based hydrogels [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]. Since ∆G enters Eq. (4.29) exponentially, substantial effects are expected following small changes in the interaction. Together with ∆G, clearly the number of nanoparticles and their size ratio with respect to the total nanoreactor size are the key quantities to tune. To save resources N should be small but large enough to warrant a high catalytic activity. The behavior of Eq. (4.29) resembles a Langmuir-binding isotherm form. The rate as a function of N initially rises linearly with a slope ε ζ exp(-β∆G) and finally saturates to the maximum rate k = k - S for large N → ∞. For a single nanoparticle, N = 1, we recover essentially the result for a yolk-shell nanoreactor k = 4πD i a exp(-β∆G), where a single nanoparticle is embedded in the center of a spherical hydrogel, apart from a slight modification of the target size, which is not a for the yolk-shell but R i , the radius of the interior hollow confinement [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF].

It is instructive to define an efficieny factor f = k/k - S between 0 and 100%, that quantifies the desired target efficiency of the nanoreactor. Solving Eq. (4.29) for N , we find

N f = e β∆G ε ζ f 1 -f (4.30)
that is, for a fixed efficiency, the number of nanoparticles needed to maintain it changes exponentially with the change of solvation free energy. As a numerical example, let us assume reasonable values of ζ = 0.2, ε = 0.01, and β∆G = -1. To obtain an efficiency of 50%, N f = 184 nanoparticle catalysts would be needed. If β∆G = -2, the number wold drop of a factor 1/e to N f = 68. Note that N f does not scale with the catalyst surface, as typical for diffusion-guided reactions, rather it decreases linearly with the catalyst size. Formula (4.30) provides a simple rule of thumb for optimizing the design and synthesis of core-shell nanoreactors. As an example, if one aims at 50% efficiency for a relatively neutral hydrogel chemical environment (∆G = 0) where the mobility of the substrate is not significantly reduced (ζ = 1), one needs to employ N f = 1/ε = R 0 /a nanoparticles. For ε = 0.01 that would be N f = 100. In the case of a polymer matrix in physical-chemical conditions leading to a reduced mobility (e.g. ζ = 0.2), one would need five times more NPs for ∆G = 0, but about the same number for β∆G ≃ -1.6. This clearly illustrates how the performance of a composite core-shell nanoreactor is non-trivially shaped by the combined action of the physical chemical properties of the hydrogel shell matrix, such as solvation free energy differences and changes in translational mobility of the substrate molecules. To exemplify some of the discussions also visually and compare the approximative Eq. (4.28) to our exact approach Eqs. (4.17), we plot in Fig. 4.3 the normalized rate versus the number of nanoparticles for two different mobility ratios ζ = 0.2 and 1.0 and core sizes γ = 0.353 (as in previous experiments described in [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]). The nanoparticle size is held fixed as ε = 0.0146 as also provided from experiments. We compare the full analytical solution Eq. (4.28) (lines) to numerical solutions of Eqs. (4.17) (symbols). First of all, we see that the analytical treatment is indeed very accurate and deviates from the exact solution by less than one percent. Comparing the different mobilities ζ within the gel, clearly all rates are higher in the more mobile case. Concerning the overall form of the curves we see the initial linear rise of the rates as predicted by the Langmuir form Eq. (4.29) and then saturating for large N values. Saturation is weak for small ∆G values but begins markedy earlier (for smaller N ) if the adsorption free energy reaches values as small as a few kBT , the thermal energy. Hence, a decisive factor in the design of optimized nanoreactors must be clearly the tuning of the reactant-hydrogel interaction towards attraction.

Analysis of the models and future perspectives

Numerous studies have been performed with the aim of characterizing the catalytic activity of active nanoreactors, and to highlight the effect of the temperature on the overall rate of the reaction. In Ref. [START_REF] Carregal-Romero | Catalysis by au@ pnipam nanocomposites: effect of the cross-linking density[END_REF][START_REF] Wu | Thermosensitive au-pnipa yolk-shell nanoparticles with tunable selectivity for catalysis[END_REF] the authors identify the interplay between the effects of the temperature on the diffusion coefficient, on the local free enthalpy ∆G and on the reactivity of the particles. The most recent analytical approach followed in [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF] models the reaction rate for a single yolk-shell nanoreactor (embedding a single nanoparticle) coupling the Debye-Smoluchowski approach with a two-state model which characterizes the swallen and collapsed state. The agreement with experimental results on yolk-shell structures accounting for single nanoparticles is good. A fundamental assumption that allows the authors in [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF] to derive an analytical expression for the rate constant with fit parameters, is that the nanoparticles do not interact with each other and can be treated independently. Our contribution allows us to describe the behavior of a composite architecture involving multiple catalysts because the equations (4.23) and (4.28) take into account explicitly the interaction among many particles. The final purpose of our study is thus to identify the values of the relevant geometrical parameters, such as size, number and configuration of the nanoparticles, that maximize the catalytic activity of a composite nanoreactor of the kind showed in Fig. 4.2. The solution worked out in the previous section allows one to compute the overall reaction rate constant of the nanoreactors for a given configuration of the N nanocatalysts (metal nanoparticles or enzymes), given the thickness of the microgel layer R 0 -R S , the diffusion coefficients in the bulk D o and in the microgel D i and the intrinsic reaction constant of the nanocatalysts k * . As in the model described in Ref. [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF], the parameters ∆G, k * , R 0 are obtained from the experimental results. However, future theoretical works should focus on describing the temperature dependence of the ligand diffusion coefficient inside the microgel and of the intrinsic rate constant, starting from first principles. Once the functional dependence of the geometrical and chemical parameters on the temperature are known, one can calculate the reaction rate constant and the approximate MOA for the catalytic system taken at a given temperature, simply by solving the system (4.17) and equation (4.23) with the corresponding parameters.

The thickness of the the microgel layer R 0 -R S varies according to changes in temperature: in the proximity of the lower critical solution temperature (LCST), which for the PNIPAM is approximately 32 • C, the shell starts to shrink due to a decreased solubility of PNIPAM in water, which causes water expulsion and compression of the porous network [START_REF] Herves | Catalysis by metallic nanoparticles in aqueous solution: model reactions[END_REF]. As a result, the shell volume is reduced and the nanoparticles get closer to each other. A functional relation between temperature and diameter of the nanoreactor can be obtained by fitting the hydrodynamic radii of the composite particle measured at different temperatures through Dynamical Light Scattering (DNLS).

In general, the intrinsic rate constant k * can be assumed to have an Arrheniuslike dependence on temperature, with pre-factor k R and activation energy E R [START_REF] Herves | Catalysis by metallic nanoparticles in aqueous solution: model reactions[END_REF]. We can thus assume, following [START_REF] Angioletti-Uberti | Theory of solvation-controlled reactions in stimuli-responsive nanoreactors[END_REF]

k * = k R exp - ∆E R k B T .
As for the value of the diffusion coefficient inside the microgel, the situation is somewhat more controversial. A plausible strategy among other seems to follow a hydrodynamic approach proposed in 1984 by Cukier [START_REF] Cukier | Diffusion of brownian spheres in semidilute polymer solutions[END_REF] for the diffusion of Brownian particles in semidilute polymer solutions with volume fraction φ

D(φ) = D 0 e -κ C φ 3/4 (4.31)
where D 0 is the diffusion coefficient at infinite dilution and κ C a microgel-dependent constant. When the temperature increases from below the critical solution temperature (LCST) and past the swollen-to-shrunk phase transition, the volume fraction of the microgel increases as it expels water molecules. For the moment, we did not include in our models the functional dependence of the parameter on the temperature, which can constitute a future development of our theory. All in all in this chapter we have demonstrated how our theory can be adapted to provide ready-to-use analytical tools to assist chemical engineers in the design of important nano-devices employed in modern applications in the chemical sciences.

Binding of small ligands to large flexible biomolecules

The interplay between competing binding sites is particularly interesting for molecules which display multi-valent activity. The coexistence of multiple binding sites can improve the efficiency of the reaction kinetics, because it increases the probability of getting in contact with a specific ligand. In the case of antibodies, a single protein can bind monovalently to two small ligands or multi-valently to larger antigens. The latter setting is used by antibodies to increase the binding affinity with large molecules and generate more stable complexes. On the other hand, we showed in the previous chapters how the competition between multiple binding sites, which mutually shield part of the flux from each other, decreases the global rate of the reaction with respect to the case of a single receptor with the same global active surface. In this chapter we aim at investigating the role of structural parameters, such as the distance and the size of active sites belonging to the same receptor, and the role of internal dynamics on the proteins' ability to bind other molecules. Dynamics is a fundamental characteristic of molecular machines, and facilitates the execution of specific biological functions [START_REF] Dennis | Molecular structure and dynamics of proteins in solution: Insights derived from high-resolution nmr approaches[END_REF][START_REF] Pa Anfinrud | Structure, dynamics, and function of proteins: New insights from time-resolved ir spectroscopy[END_REF]. As underlined by many experimental data, proteins can not be considered as rigid and static object because they exhibit a variable degree of structural flexibility. Flexible units might act as dynamical gates that govern the accessibility of specific sites and indirectly control the cascade of reactions triggered by a binding event [START_REF] Dennis | Molecular structure and dynamics of proteins in solution: Insights derived from high-resolution nmr approaches[END_REF][START_REF] Pa Anfinrud | Structure, dynamics, and function of proteins: New insights from time-resolved ir spectroscopy[END_REF].

A paradigmatic example of the proteins' complexity is given by the typical antigen-antibody reaction. Antibodies are large, extremely flexible molecules with multiple binding sites, whose internal dynamics is certainly key to their great ability to bind antigens of all sizes, from small hormones to giant viruses. Unarguably, they can be counted among the most important molecular machines for the functioning of life. Like other large biomolecular assemblies, they are increasingly being exploited in modern nanobiotechnology [127] and biomedical [START_REF] Justin | The history of monoclonal antibody development -progress, remaining challenges and future innovations[END_REF] applications. Hence, a satisfactory comprehension of antibodies' structure and functional dynamics poses great challenges at the interface among biology, physics, medicine and pharmaceutical industrial research. We will base our analysis on the structure of a particular class of antibodies, the immunoglobulin G (IgG), which is the isotype that provides the majority of antibody-based immunity against invading pathogens. Antibodies are large Y-shape proteins whose ultimate function is to bind hostile organisms, such as viruses and bacteria. The reaction between an antibody and an antigen occurs through one (or both) of the two active sites placed on the IgG surface, which bind to the antigen's epitope(s). As we will precise in the following section, the different domains of an antibody are connected together in a non-rigid fashion, allowing the molecule to adopt extremely variable configurations. The question we investigate is how the internal dynamic of the antibody influences the efficiency of the binding reaction.

In this section, we first build a shape-based coarse-grained model of an IgG molecule and then we employ it in a molecular dynamic simulation to generate plausible 3D conformations of the molecule in agreement with experimental data (single-molecule Cryo-ET ). Then, we apply our theoretical model, introduced in the previous chapter, to compute the binding rate constants of small antigens to the IgG which is assumed frozen in each one of the obtained three-dimensional conformation. Therefore, we implicitly assume that the ligand (antigen) molecules are small enough, so that a local equilibrium can be reached faster than the typical rearrangement times of internal IgG conformations. Finally we compare the results corresponding to different configurations and we identify the meaningful structural parameters for the modulation of the rate. Moreover, we map the problem of the binding between antibody and antigen onto a simple effective model, which can be described with very a reduced set of parameters. Our findings pave the way for further investigation of the subtle connection between the dynamics and the function of large, flexible multi-valent molecular machines.

The IgG and its representation

Immunoglobulins are large glycoproteins of approximately 150kDa made of four polypeptidic chains: two light chains (of ∼ 25 kDa) and two heavy chains (of ∼ 50 kDa) which give the IgG a characteristic Y shape, which is shown in fig 5 .1. In the macroscopic structure of the immunoglobuline it is possible to recognize three roughly ellipsoidal subunits of length about 6 nm: two fragment binding arms (Fabs), and a stem (Fc), connected by a flexible hinge [START_REF] Harris | The three-dimensional structure of an intact monoclonal antibody for canine lymphoma[END_REF][START_REF] Harris | Comparison of the conformations of two intact monoclonal antibodies with hinges*[END_REF][START_REF] Ollmann Saphire | Contrasting IgG structures reveal extreme asymmetry and flexibility[END_REF].

The tips of the Fab domain host hypervariable regions [START_REF] Elgert | Immunology: Understanding the Immune System[END_REF] (see Fig. 5.2). These are often referred to as the active sites or paratopes, as these are the portions of the structure where antigens are bound (at their epitopes). The Fc stem is recognized at its lower end by the complement system [START_REF] Van Lookeren Campagne | Macrophage complement receptors and pathogen clearance[END_REF] and by phagocytic cells [START_REF] Lynda | Phagocytosis: elegant complexity[END_REF] in the early steps of an immune response. These three lobes are hinged together by a flexible stretch of polypeptide chain, which allows them to vary considerably their mutual positions.

Atomic-force microscopy [START_REF] Kienberger | Following single antibody binding to purple membranes in real time[END_REF] and cryo-microscopy [START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF] measurements have revealed that Fab-Fab and Fab-Fc angles are virtually limited only by steric clashes, with measured values ranging from 15 o to 128 o (Fab-Fc angle) and from about 20 o to 180 o (Fab-Fab angle) [START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF].

The intrinsic flexibility of IgG molecules reflects their ability to bind antigens of different sizes, from small molecules such as hormones to large viruses [START_REF] Dennis R Burton | Antibody and virus: binding and neutralization[END_REF][START_REF] I S Mian | Structure, function and properties of antibody binding sites[END_REF]. Moreover, high flexibility is also key to double-Fab (bivalent) binding to large viruses [START_REF] Wang | Neutralization efficiency is greatly enhanced by bivalent binding of an antibody to epitopes in the v4 region and the membrane-proximal external region within one trimer of human immunodeficiency virus type 1 glycoproteins[END_REF], a process quantified by the so-called binding avidity [START_REF] Ns Greenspan | Affinity, complementarity, cooperativity, and specificity in antibody recognition[END_REF][START_REF] Mathai Mammen | Polyvalent interactions in biological systems: Implications for design and use of multivalent ligands and inhibitors[END_REF], as opposed to the single-binding affinity. Bivalent binding increases the overall strength of the immune response and also allows for IgG-mediated virion aggregation [START_REF] Sam | Valency of antibody binding to enveloped virus particles as determined by surface plasmon resonance[END_REF]. Re-cently, it was demonstrated that the intrinsic flexibility of antibodies can be exploited to have them literally walking on antigen-covered surfaces with specific lattice-like arrangements of haptens with lattice spacing matching the IgG stride [START_REF] Preiner | IgGs are made for walking on bacterial and viral surfaces[END_REF].

In this chapter we focus on the following question: given their great flexibility, we want to assess whether IgGs are more effective in binding a co-diffusing antigen when adopting certain specific conformations. More generally, it would be extremely interesting to establish a quantitative link between the large-scale dynamics of antibodies and their binding efficiency. Here we concentrate on this problem in the case of small antigens, where the dynamics of substrate and IgG molecules are characterized by widely separated time scales. More precisely, we elaborate a model to assess the role of the IgG conformation on the diffusion-limited reaction rate describing the formation of an encounter complex between a diffusing antigen and a stationary IgG molecule frozen in a given conformation.

To obtain an effective model of the antibody, we preliminary consider its static atomistic representation, as obtained via X-ray crystallography. More specifically, we refer to the intact map of the murine IgG 1IGT as made available by the RCSB protein data bank, (transparent structure in figure 5.3). We built a coarse grained model of the immunoglobulin made of N = 96 identical spheres of radius R, automatically docked into the crystal structure by an ad hoc fitting routine made available by the computer software VMD [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF]. The diameter of the beads composing the skeleton of the antibody is equal to the smallest distance between pairs of spheres, multiplied by 0.95, a scaling factor which prevents unphysical merging of contiguous units. Based on this criterion, we obtain R=4.4 Å. The radius of the active sphere, localized on the outer edges of the Fab is assigned different values (in the following the size of the paratopes will be indicated in the captions and figures' legends). From crystallographic data, the active surface for the antigen-antibody complex seems to vary with the antigens size. It scales linearly for very small antigens, but then saturates to an almost constant value which never exceeds 1200 Å. Our choice for the number of beads in the coarse grained representation is a good trade-off between the constraints of keeping the structure light enough (fewer spheres) and reproducing faithfully enough the three-dimensional shape of the domains (more spheres). As we will state in the following, the main conclusion of our study holds in general, irrespectively at the specific level of coarse graining imposed in the model.

An illustrative sketch of the coarse grained model is depicted in Fig. 5.3. The active sites, situated on the top of the Fab arms, are represented by two absorbing sinks of appropriate size, while the spheres which define the bulk of the structure are assumed to be totally reflecting. The above structure can be readily turned into a full dynamical model to reproduce the internal dynamic of the protein. We will describe hereafter the molecular dynamic simulation that we developed to obtain a set of independent configuration of a single IgG molecule. Following the well known elastic-network strategy [RFF] we map our IgG into a network of beads and springs. This is achieved by stretching a Hookean spring of stiffness κ between any two beads that are separated by less than a specific cutoff length R c in the equilibrium structure (by definition, the X-ray structure depicted as the transparent structure in Fig. 5.3). The constant κ is chosen stiff enough so that the shape of each domain is preserved while they fluctuate freely about one another. In order to prevent steric clashes among the three lobes, we also introduce a repulsive potential energy between any [START_REF] Arkhipov | Stability and Dynamics of Virus Capsids Described by Coarse-Grained Modeling[END_REF][START_REF] Arkhipov | Four-scale description of membrane sculpting by BAR domains[END_REF] implemented in the VMD package [START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF]. The black spheres positioned on the edges of the Fab represent the active sites. In this schematic illustration, the beads occupy the positions that are assigned through docking. The radius of the beads shown in the figure corresponds approximately to the one used in the simulations (0.44 nm). The atomic representation of the antibody is shown as a transparent structure (PDB code: 1IGT) two beads. The total potential energy reads:

V = κ 2 i>j c ij (r ij -R ij ) 2 + v 0 i>j ℓ r ij 12 .
(5.1)

Here the position of the i-th bead as a function of time and in the equilibrium structure are denoted by r i (t) and R i , respectively. Accordingly, R ij = |R j -R i | and r ij (t) = |r j (t) -r i (t)| are the equilibrium and displaced inter-bead distances, respectively. The matrix

c ij = 1 if R ij < R c 0 otherwise
specifies all the interacting pairs, and is known as the connectivity matrix. In the following we take v 0 = 5 kcal/mol, ℓ = 10 Å, κ = 50 kcal/mol/ Å2 and R c = 22 Å. A set of different configurations was then obtained by sampling a constant-energy trajectory with initial conditions given by the crystallographic structure and initial velocities assigned randomly to the centers of the beads. In order to reconstruct the trajectory, we integrated Newtons equation numerically

mr i = -∇ i V
through a position-extended Forest-Ruth like (PEFRL) symplectic algorithm [START_REF] Ip Omelyan | Optimized forest-ruth-and suzuki-like algorithms for integration of motion in many-body systems[END_REF] with a time step dt = 2.8 fsec. The mass m of the effective beads was fixed at m = M IgG /96 ≈ 1.67 kDa. The trajectory was sampled every 0.56 psec to obtain many sets of coordinates, providing the skeletons of different conformers. The snapshots in figure 5.4 displays three samples of our trajectory.

In order to validate the results of our Coarse-grained molecular dynamics simulations, we compared the statistics of the inter-domain angles with the results of Cryo Electron Tomography (Cryo-ET) measurements, taken from Ref [START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF]. In [START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF] an ensemble of Cryo Electron Tomography experiments (see fig 5.5) is analyzed to reconstruct the distribution of configurations visited by single IgG molecules in solution. To interpret the reconstructed three-dimensional portraits, a simplified representation of the immunoglobulin was put forward: the IgG was sketched as composed of three rigid rods freely jointed together in a common point, the hinge. By fitting such mechanical model to the experimental structures, the Fab-Fc angle and the Fab-Fab angle for all the reconstructed molecules were registered and their frequency of occurrence determined. The obtained cumulative distributions for φ and ψ are depicted in figure 5.6 and compared to the results of our simulations. The Fab-Fab and Fab-Fc angles were obtained for each conformer resulting from our dynamical model by first constructing the inertial ellipsoids of the three domains. The three angles (ψ and the equivalent φ 1 , φ 2 ) were then computed from the scalar products of the vectors describing the major axes of the three ellipsoids (pointing outward from the hinge). The surprising conclusion of the comparison is that a simple shape-based coarse-grained model that is able to reproduce the correct shape of the antibody suffices to reproduce the experiments in solution. We stress that our simulations are deterministic and carried out in vacuum: the IgG's domains The image is reproduced from Ref. [START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF].

fluctuate about the hinge (the center of mass is at rest) only subject to mutual collisions, while preserving their shape due to the stiff springs stretched among the beads. The reason behind this somewhat surprising finding is that we are looking at equilibrium properties. It is then manifestly redundant what the collisions with the solvent add to the random collisions among the domains -they would only generate an equivalent noise spectrum that would not change the equilibrium statistics. This also shows that hydrodynamic effects seem not to affect to an appreciable extent the large-scale structural fluctuations of IgGs, which appear mainly controlled by the excluded-volume effects related to the shape of the mutually hinged domains.

Mathematical model for the binding reaction

The collection of configurations obtained with the dynamical coarse grained model represents the starting point for our analysis through the method introduced in chapter 3. The chemical reaction between an antibody and an antigen can be represented as a two-step process, according to the scheme

[IgG] + [A] k D ---⇀ ↽ --- k -D [C] * k * -→ [C]
which describes first the formation of the antigen-antibody complex, and eventually the stabilization of the antigen-antibody complex.

If we restrict our scenario to small antigens, which are diffusing sufficiently fast so to see the antibody virtually frozen in one of the allowed configuration, and we assume that antibodies are sufficiently diluted when compared to the concentration of the other reactants, we can map the chemical reaction on a boundary problem of the kind (3.6). More precisely, assuming a separation between the characteristic time scale for the diffusion of the antigens and for the time scale associated with the IgG's internal dynamics, we define a set of boundary problems, one for each configuration obtained through the MD simulation. For a given configuration, the IgG is represented as a collection of disconnected spherical boundaries Ω α (α = 1, ..., N ) placed at given locations in space (see Fig. 5.7). The volume external to the Ω α spheres is imagined to be filled by a continuous concentration u of diffusing antigens, which are moving with an effective diffusivity D = D A + D IgG and are characterized by a bulk concentration c 0 (far from the antibody). The IgG configuration is hidden in the boundary conditions. The outer spheres at the Fab tips of radius R a , meant to represent the paratopes, are considered as perfectly absorbing boundaries. The remaining N -2 spheres with equal radius R are taken to be totally reflecting. With all these assumptions, the steady-state (t → ∞) rate constant k corresponding to a given conformation of the IgG can be calculated analytically by computing the total flux to the spheres which represent the paratopes (Ω 1 and Ω 2 hereafter). The stationary problem that we have to solve for the normalized density u = c/c 0 of the antigens reads The solid line stands for the results of our coarse-grained simulations, while symbols refer to the Cryo-ET experiments reported in Ref. [START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF]. The inset in panel (a) provides a schematic view of the mechanical model employed in the simulation, together with a definition of the relevant angles ψ and φ.

               ▽ 2 u =
The system (5.2) is solved by employing the method described in the previous sections. The solution can be written formally here as combination of only irregular harmonics, since we are working in an unbounded domain:

u = 1 + N α=1 u α = 1 + N α=1 ∞ n=0 n m=-n B α mn u - mn (r α ) = 1 + N α=1 ∞ n=0 n m=-n B α mn r α R α -n-1 P m n (cosθ α )e imφα .
(5.

3)

The total antigen flux across the paratopes surfaces, normalized to twice the reaction rate of an isolated spherical sink with the same radius R a as the active paratope spheres, reads

k 2k S = - (B (1)
00 + B

(2) 00 ) 2 .

(5.4)

In Fig. 5.8, we plot the normalized rate constant k as a function of the distance between the two binding sites (Fab tips), after ordering the different configurations in terms of the relative distance d between the two paratopes. Several comments are in order after examining the results. As a first observation, the global rate of the antibody-antigens reaction is always smaller than twice the rate of an isolated paratope sink. This reduction stems from the competition of the paratopes for the target antigen molecule and from non trivial screening effects exerted by the reflecting IgG body. The conclusion that we can draw is that the interaction among the sinks is always anticooperative, and the rate to a selected paratope increases as the distance between the sinks grows. The configurations which make the binding reaction more efficient are thus the ones associated with the maximum value of the Fab-Fab angle. We will see in the following how the configurations which result in a higher binding efficiency are also the most visited ones (see again the second panel in figure 5.6).

The second remark is that the points displayed in Figure 5.8 fall on a smooth curve. The absence of appreciable scattering around the average profile seems to suggest that the relative distance between the paratopes is indeed the only meaningful parameter that needs to be retained when aiming at a quantitative characterization of the binding performance of a dynamical IgG. The reflecting beads of the IgG body shield out part of the incoming flux of antigens. However, such screening seems not to depend on how the inert spheres are arranged around the paratopes, so long as the paratopes are at a fixed distance.

Moreover, figure 5.8 shows that the results can be interpreted by mapping the antibody to an effective and rather simple model: the IgG can be replaced by a dumbbell made of two spherical sinks of radius R a placed at a distance d. The red curves in the figure represent in fact the analytical expression of the reaction rate in the monopole approximation for each of two absorbing identical spheres (3.53), multiplied by a constant factor f a < 1 which accounts for the excluded volume effect of the reflecting beads:

k 2k S = f a 1 1 + R a /d , (5.5)
where k S = 4πDR a is the rate constant of an isolated paratope (see caption of Fig. 5.8). The factor f a is the only remnant of the IgG body (the other N -2 reflecting spheres), whose action on the two sinks only causes a rather modest 6 ÷ 7 % reduction with respect to an isolated paratope dumbbell.

To shed further light on the above findings, it is instructive to look at the distributions of the computed reaction rate constants. These are shown in Fig. 5.9. One can appreciate how the flexibility of the IgG impacts the statistics of the rate constants. The profiles appear negatively skewed, which implies that the higher weight is associated with larger values of the rate constant. Interestingly, the histograms can be fitted to Gumbel profiles, which suggests that the rate constants computed from the sampled configurations may correspond to near-maximum sampled values. Moreover, the curves shift towards the right, namely to higher values of the normalized rate constants, when R a is reduced. This is a consequence of the reduced diffusion interaction between the two active sites. In fact, the larger the active regions, the greater the rate constant for an isolated paratope, but the smaller the combined rate constants for the two-paratope system [START_REF] Traytak | The diffusive interaction in diffusion-limited reactions: the steadystate case[END_REF][START_REF] Piazza | Anticooperativity in diffusion-controlled reactions with pairs of anisotropic domains: A model for the antigen-antibody encounter[END_REF][START_REF] Bongini | A dynamical study of antibody-antigen encounter reactions[END_REF]. Therefore, we conclude that internal flexibility forces the antibody to visit preferentially those configurations that are associated with higher probability of forming an encounter complex. Summarizing, our results suggest that large, flexible molecules may have been designed by evolution to exploit their flexibility to a maximum degree in terms R a = 0.44 nm R a = 0.54 nm of the ability of binding small antigens. Of course, this is only part of the story. IgGs bind also bivalently to different epitopes on the same surface, such as a virus capsid. Thus, the mutual flexibility and shape of the three lobes must also have been shaped by the evolutionary pressure exerted by requiring such binding events to be optimized.

The results illustrated so far hold in general, and do not depend on the level of coarse-graining imposed in building the effective bead-based model of IgGs. To prove the robustness of this strategy, we want to show how the same analysis can be repeated for different choices of the number of beads that make up the coarsegrained IgG molecule, yielding similar conclusions. However, one could in principle speculate that a sphere-based model will inevitably yield a structure with plenty of holes, unphysical consequences of the coarse-graining procedure. The antigens, assumed point-like as their size is incorporated in the size of the paratope beads, R a , could diffuse through the structure and this could be the reason why the whole IgGs behaves quite just like a nearly perfect dumbbell. In order to investigate further into this matter, we studied a hyper-simplified model where the IgG is replaced by three contacting spheres, two representing the Fab arms and one modelling the Fc stem. The active paratopes were still modelled as two additional absorbing spheres at the outer ends of the Fab spheres (see cartoons in Fig. 5.10). Fig. 5.10 illustrates our results for this ultra-coarse-grained model. Overall, the picture traced in the previous paragraph is confirmed, which shows that the partially transparent nature of the multi-bead model does not introduce artifacts in the computation of the rate constant and can be safely adopted as it allows one to sample the large-scale configuration space in agreement with the experiments. The rate constant in the three-bead model can still be described by a modified (rescaled) dumbbell: the largest rate constants are invariably associated with the configurations where the Fab arms are stretched away from each other to a maximum, while the smallest rate constants correspond unfailingly to the situation where the active sites are close to each other (see cartoons in Fig. 5.10). The rescaling factor to be used in the dumbbell law is close to one for very large and very small sizes of the paratopes. In the former case, the reason is that the active spheres are strong enough not to feel the screening of the three-bead reflecting IgG. In the latter case, the paratopes are so small that their separation is always several times their size. At intermediate sizes of the active sites the reduction is largest.

To summarize, we have developed a coarse-grained dynamical model of the IgG which has been exploited to study the interplay between the inherent dynamics of the molecule and its specific functioning. The first result of our analysis is that simulations of the coarse grained model of the IgG in vacuum reproduce with good accuracy the statistic of the main structural parameters, as extracted from single molecule Cryo-ET experiments. We then set down to quantify the reaction rate between the dynamical coarse-grained IgG and a sea of, freely diffusing, small antigens, inside a given volume, subject to open boundary conditions. The study exploits the solution of the Laplace equation in spherical coordinates as it is carried out in the previous chapter. Working in this framework, we showed how our method, together with a coarse grained strategy, allows to extend the simpler Smoluchowski model to the case of complex shaped large biomolecules. The spherical representation is in fact too simplistic for many biological molecules: we showed how we can compute the rate of a reaction for an arbitrary choice of the shape of the receptor, by reproducing its configuration with assemblies of spheres. The obtained results can then be rationalized by resorting to an elementary interpretative model, which assumes dealing with just two sinks. The fact that the behavior of complex three dimensional objects can be nicely explained by simplified dynamical models, constitutes another interesting conclusion of this study. Moreover, we underlined how the internal flexibility forces the antibody to visit preferentially that subset of configurations which is associated with the larger values of the predicted reaction rate. We interpreted this correlation as a sign of the efficiency of the protein, which nature could have engineered so as to minimize the negative interference between its two active sites. The internal flexibility keeps apart, on average, the two paratopes and thus enhances the performance of their cooperative binding. Of course, our kinetics results are theoretical, and it would be very interesting to compare our predictions to direct experimental measurement on the kinetic of the IgG. Unfortunately, the effect of the internal dynamics on the kinetics of the reaction is hard to detect experimentally, since conventional measurements return global estimates and do not allow to disentangle the information on specific IgGs configurations. One solution could be to compare measurements performed with the same antigens on wild-type antibodies and on specifically designed molecules, whose conformational ensembles are confined to selected portions of their phase space. This might be the case of antibodies designed to partially reduce or, even delete the hinge domain, thus inducing a substantial loss of flexibility [START_REF] Guddat | Three-dimensional structure of a human immunoglobulin with a hinge deletion[END_REF]. An alternative route could be to compare the kinetics of (i) individual antigen-binding fragments (Fabs) and (ii) assemblies of two, three or more recombined such fragments with the kinetics of intact antibodies [START_REF] Todorovska | Design and application of diabodies, triabodies and tetrabodies for cancer targeting[END_REF]. Another interesting and viable route to test the predictions of our theory would be to analyze antibodies from camelids, such as dromedary or llamas [START_REF] Desmyter | Camelid nanobodies: killing two birds with one stone[END_REF]152]. About half of the antibodies of these animals lack light chains and only feature two heavy chains with three IgG domains each [152], with the missing IgG domain replaced by a flexible linker. The IgG domain carrying the complementarity determining regions (CDR) display two point mutations that make their surface more hydrophilic (where normal IgG feature the interface between the terminal light-chain and heavy-chain CDR-carrying domains). This has prompted researchers to isolate these special IgG domains, which are now known under the name of nanobodies and hold great promise for the biotechnology industry [START_REF] Desmyter | Camelid nanobodies: killing two birds with one stone[END_REF]. Interestingly, nanobodies show the same high affinities of full antibodies and can be used to construct more complex molecules, such as by joining two of them via a long and flexible hinge. It would be interesting to compare kinetics measurements performed on a host of different antigen-binding systems, such as the one described above, to validate our theoretical approach.

As a final remark, it is worth emphasizing that the techniques here developed are general and could be applied to study complex formations involving a wide host of different biomolecules of various sizes and flexibility.

Chapter 6

Gated reactions

The fundamental assumption underlying the model described in the previous chapters is that the system relaxation time for the diffusive flux of particles (the ligands) t D is small enough to neglect time-dependent effects. Hence, for the applications presented, we were able to consider a local equilibrium of the ligand concentration around the reactive boundaries and hence restrict ourselves to the stationary diffusion equation. Often though, fluctuations in the conformation and reactivity of the ligand or the receptor play an important role in the extent of the binding reaction. The classical example is the binding of oxygen and myoglobin which could not happen if the protein were rigid, since the binding site can be occluded as a result of the protein's conformational dynamics. The same kinetics is observed in the binding of acetylcholinesterase and acetylcholine, where the active site is placed at the end of a bottleneck which opens and closes responding to thermal fluctuations [START_REF] Andrew | Gated diffusion-controlled reactions[END_REF]. Processes like these can not be described by stationary models, because the temporal fluctuations occur on a time scale comparable, or even smaller than the characteristic time for the diffusion of the molecules.

In this chapter we will consider problems where the boundary conditions are time dependent. We will use a different analytical approach, since the harmonicbased technique used up to now can not be exploited directly for solving the time dependent diffusion equation. We will rely on Duhamel's theorem, a powerful theory introduced by the French mathematician Jean-Marie Duhamel. Duhamel provided a convenient approach for developing a solution to diffusion equations with time-dependent boundary conditions and equations with a time dependent energy generation term [START_REF] Necati | Heat conduction[END_REF]. The idea behind his method is to exploit the solution of the corresponding problem with time-independent boundary conditions, to build the solution of the original problem. This procedure is described in Duhamel's principal work on heat conduction. We will adopt here the same representation of the ligand /receptor interaction introduced in chapter 3, where the hypotheses of high dilution, little diffusivity and small concentration of one reactant allowed us to represent it as a fixed sphere in the reference frame and to define a boundary problem for the density of the diffusing species. For the sake of simplicity, we will restrict here to the Smoluchowski setting, namely the approximation in which the encounter probability of two receptors is neglected and we can consider a single particle Ω of the species A within an unbounded domain.

To provide a description of a gated reaction we have to take into account the time dependent equivalent of system (3.3), complemented by time-dependent boundary conditions on the encounter surface:

         ∂ρ ∂t = D∇ 2 ρ ρ(r, 0) = ρ 0 (r) r ∂ρ ∂r -h(t)ρ = 0 on ∂A. (6.1)
Here ρ(x, t) is the density of a population which is diffusing in the open space Ω = R 3 \ A external to a spherical boundary of radius R centered in the origin of the reference system. The radius R of the sphere is an effective dimension which represents the encounter distance between the reactants, and, similarly, the coefficient D embodies an effective diffusivity, as introduced in chapter 3. The Laplacian operator, given the radial symmetry of the problem, reduces to

∇ 2 ρ = 1 r 2 ∂ ∂r (r 2 ∂ρ ∂r )
.

The function h(t) in the boundary condition gauges the fluctuations in the reactivity of the boundary and it accounts for the typical time variation of the absorbing properties of the receptor. As we will see after deriving the solution, the function h(t) will have to match some regularity conditions to make the analytical calculations possible.

The previous system can be reduced to a one-dimensional problem for u = rρ by exploiting the radial symmetry of the domain :

         ∂u ∂t = D ∂ 2 u ∂r 2 in (R, ∞) u(r, 0) = rρ 0 (r) in (R, ∞) ∂u ∂r = g(t) for r = R, (6.2)
where R is the effective radius of the receptor and the function g(t) allows us to express the boundary conditions as Neumann boundary conditions:

g(t) = 1 R h(t) + 1 u(R, t) = h(t) R u(R, t). (6.3)
This trick is necessary to employ Duhamel's theorem. To understand better the philosophy of Duhamel's approach, we will briefly describe the general statement of the theorem for a diffusion equation subject to time-dependent boundary conditions

         ∂T (x, t) ∂t = D∇ 2 T (x, t) in Ω T (x, 0) = F (x) in Ω ∂T ∂n i + h i T = λ i (x, t) on ∂Ω i . (6.4) 123
by means of the solution of the auxiliary problem:

           ∂ϕ(x, t, τ ) ∂t = D∇ 2 ϕ(x, t, τ ) in Ω ϕ(x, 0, τ ) = F (x) in Ω ∂ϕ(x, t, τ ) ∂n i + h i ϕ(x, t, τ ) = λ i (x, τ ) on ∂Ω i . (6.5)
Here t is the time, while τ is an auxiliary parameter which allows one to make the functions λ i (x, τ ) not dependent on time. The relation between F and ϕ is given by

T (x, t) = F (x) + t τ =0
∂ϕ(x, t -τ, τ ) ∂t dτ. (6.6)

Once the system (6.5) is solved for ϕ(x, t, τ ), it is thus possible to obtain the solution of (6.4) using equation (6.6).

In principle, this approach can not be used directly to solve the gated boundary problem (6.2), because the time dependency of our problem is hidden inside the intrinsic rate constant h(t) of the radiative boundary condition. This is the reason why in equation ( 6.2) we have introduced the function g(t) and we have formally written the boundary condition as Neumann's. This strategy will enable us to solve the gated boundary value problem, even if the analytical expression of the solution will be given implicitly.

Let then us focus on applying the Duhamel theory to the problem (6.2). After introducing the parameter τ , the auxiliary problem reads:

         ∂φ ∂t (r, t, τ ) = D ∂ 2 ∂r 2 φ(r, t, τ ) in (R, ∞) φ(r, 0, τ ) = rρ 0 (r) in (R, ∞) ∂φ ∂r = g(τ ) for r = R. (6.7)
We remind that here g(τ ) is not a function of time. In order to solve system (6.7) we decompose the solution as φ(r, t, τ ) = φ 1 (r, τ ) + φ 2 (r, t, τ ), where φ 1 obeys to the steady state equation:

     ∂ 2 φ 1 (r, τ ) ∂r 2 = 0 in (R, ∞) ∂φ 1 ∂r = g(τ ) for r = R. (6.8)
while φ 2 is the solution of:

         ∂φ 2 ∂t (r, t, τ ) = D ∂ 2 φ 2 (r, t, τ ) ∂r 2 in (R, ∞) φ 2 (r, 0, τ ) = rρ 0 (r) -φ 1 (r, τ ) in (R, ∞) ∂φ 2 ∂r = 0 for r = R. . ( 6
.9)

Thus we have the following solution for the auxiliary problem (6.7):

φ(r, t, τ ) =φ 1 (r, τ ) + φ 2 (r, t, τ ) =g(τ )r + 1 √ 4πDt ∞ x=R (xρ 0 (x) -g(τ )x) exp - (r -x) 2 4Dt + exp - (r -2R + x) 2 4Dt dx (6.10)
By using the relation (6.6) between the solution of the auxiliary problem and the solution of the time dependent problem, we have:

u(r, t) =rρ 0 (r) + t 0 ∂ ∂t 1 4πD(t -τ ) ∞ x=R xρ 0 (x) - h(τ ) R u(R, τ )x exp - (r -x) 2 4D(t -τ ) + exp - (r -2R + x) 2 4D(t -τ ) dx dτ =rρ 0 (r) - t 0 ∂ ∂τ 1 4πD(t -τ ) ∞ x=R xρ 0 (x) exp - (r -x) 2 4D(t -τ ) + exp - (r -2R + x) 2 4D(t -τ ) dx dτ - t 0 h(τ ) R u(R, τ ) ∂ ∂t 1 4πD(t -τ ) ∞ x=R x exp - (r -x) 2 4D(t -τ ) + exp - (r -2R + x) 2 4D(t -τ ) dx dτ = 1 √ 4πDt ∞ x=R xρ 0 (x) exp - (r -x) 2 4Dt + exp - (r -2R + x) 2 4Dt dx - t 0 h(τ ) R u(R, τ ) ∂ ∂t 4D π (t -τ ) exp - (r -R) 2 4D(t -τ ) +(r -R)erf (r -R) 4D(t -τ ) dτ.
(6.11)

The solution of (6.2) obtained here is implicit in the sense that the expression of u(r, t) involves u(R, t) and thus it can not be computed directly. Remarkably, we notice that, for our purpose of computing the reaction rate

k = D δΩ ∂ρ ∂r dσ,
we only need to know the value of the solution restricted to r = R. In fact, after substituting the boundary condition in (6.2) in the definition of the normalized rate constant one gets: .12) where k S = 4πDRρ B is the stationary Smoluchowski rate of an absorbing ungated sphere and ρ B is the bulk density. Thus, to get the reaction rate, we need to evaluate the implicit expression (6.11) for r = R. This leads to:

k(t) k S = ∂ρ(R, t) ∂r = h(t) ρ(R, t) ρ B . ( 6
u(R, t) = l(t) - t 0 K(t, τ )u(R, τ )dτ (6.13)
where

K(t, τ ) = h(τ ) D πR 2 (t -τ ) l(t) = 1 √ πDt ∞ R xρ 0 (x) exp - (x -R) 2 4Dt dx. (6.14)
Expression (6.13) is a linear Volterra integral equation of the second kind. Solutions of Volterra equations of the second kind can be found by Laplace transformation only for convolution kernels K(t, τ ) = K(tτ ). Since in our case the kernel involves a convolution function multiplied by a function of the variable t only, h(t) = (h(t) + 1), the equation is solvable analytically only for a choice of h(t) = const.

It is instructive to check (6.13) with h(t) = h, which amounts to consider a timeindependent radiative boundary condition. It is not difficult, by using the Laplace transform, to see that this choice inserted in (6.13) gives in the stationary limit

k k S = h 1 + h , as expected.
Even though this choice allows us to compute the solution analytically, it is not of great interest for our purposes, because the gated behavior of the boundary is lost in the choice h(t) = const. Among other strategies for solving the Volterra integral equation, simple algorithms based on quadrature schemes can be used to solve iteratively equation (6.13) when the kernel K(t, τ ) and function l(t) are continuous for all 0 ≤ τ ≤ t.

The main difficulty here is the fact that the kernel in (6.14) shows a weak discontinuity for t = τ . The Volterra integral equation resulting from discontinuous but integrable kernels involving functions of the kind (ts) α f (t, s) 0 < α < 1 is called Weakly Singular Volterra Integral Equation. Specific methods have been proposed by several authors for equations of this kind. Here we make use of the algorithm developed using a fractional backward differentiation formulae (BDF) method described in Ref. [START_REF] Lubich | Fractional linear multistep methods for abel-volterra integral equations of the second kind[END_REF][START_REF] Hairer | Fast numerical solution of weakly singular volterra integral equations[END_REF][START_REF] Baker | Fft techniques in the numerical solution of convolution equations[END_REF]. This method allows one to solve a Weakly Singular Volterra-Abel Integral Equation of the second kind exhibiting the following structure: y(t) = s(t) + 1 π t 0 k(tτ ) √ tτ g(τ, y(τ ))dτ, 0 ≤ t ≤ T. (6.15) This method requires that the functions involved are sufficiently smooth. This scheme applies directly to equation (6.13) with k(tτ ) = -D and g(τ, y(τ )) = h(τ )y(τ ). Once we have defined the specific boundary value problem, by choosing a gating function h(t) for the boundary condition and the initial distribution of the density f (r) in (6.1), then the solution u(R) is computed from (6.13) and the normalized time dependent rate is:

k(t) k S = 1 Rρ B h(t)u(R, t).
Together with the time dependent rate, we will compute for every choice of the parameters the mean rate constant k which is the time-dependent rate constant averaged over a typical integration time T , that is selected according to the properties of the gating function:

k = 1 T T 0 k(t)dt.

Gated reactions in the literature

The internal dynamic of a molecule can have a strong effect on the extent of a reaction with another agent if the time scale associated with the conformational or reactivity fluctuation is comparable with the characteristic time of diffusion. In our previous analysis of the antibody/antigen interaction (see section 5) we worked under the hypothesis that the antigens diffuse sufficiently fast with respect to the typical time of large-scale conformational rearrangements of an IgG molecule. For a broad choice of the parameters, assuming that a ligand sees the receptor frozen in a given conformation (both geometrical and/or relative to the chemical properties) is not a realistic assumption. In those cases, a time dependent description which takes into account the role of the reactivity fluctuations has to be adopted. A quantitative analysis of the kinetic effects of the gated binding processes was first carried out by McCamon and Northrup [START_REF] Mccammon | Gated binding of ligands to proteins[END_REF], who identified two limiting situations for proteins (such as lysozyme) whose active sites are either accessible or not due to internal motions of the molecule. The first limit is the fast gating limit, for which opening and closing of the gate is much faster than the time required by a ligand to escape from the neighborhood of the gate. The second is the slow gating limit, in which the protein stays for a long time in the open and in the closed state compared to the rate of escape of the ligand. They restricted their study to sinks which are totally absorbing when in the open state and which require negligible time to open or close.

As an application of this analysis, we recall the recent study in [START_REF] Wade | Simulation of enzymesubstrate encounter with gated active sites[END_REF] in which Wade et al. use a dynamical simulation to analyze the effect of fluctuations in reactivity of ligand-receptor binding. They focus on the kinetic analysis of chicken muscle triosephosphate isomerase, an enzyme which binds to glyceraldehyde 3-phosphate through an active site which shows flexible loops that act as an oscillating gate.

The opening/closing of the gate occurs with an average period of 1 ns while the characteristic time of diffusion for the ligand is estimated as 16 ns. Given these estimates, the kinetics falls within the fast gating regime. Both simulations and experiments show that the gating has a very little effect on the rate of the enzyme, which is diffusion controlled. In all the mentioned studies, as well as in the majority of the works on this topic, gated reactions have been studied mainly stochastically, by taking the fluctuations of the gated active site as Markov processes. In this perspective the reaction can be described by the following scheme:

[open]

a -⇀ ↽ - b [closed].
In Ref. [START_REF] Szabo | Stochastically gated diffusion-influenced reactions[END_REF] it is shown that in the regimes of fast gating, where both the rates a and b are very large, the system behaves as if the gating did not exist, so that the overall rate of the binding reaction is approximately the same of a constantly reactive receptor. When the opening and closing of the gate is slow compared to the characteristic time for diffusion, so that (a + b) -1 ≫ τ D , the gated rate constant is

k = b a + b k ∞ (6.16)
that is the steady state reaction rate of a totally absorbing sink multiplied by the probability of finding the gate open. Nevertheless there exist some studies which employ deterministic gating functions, as in Ref. [START_REF] S H Northrup | Rate theory for gated diffusioninfluenced ligand binding to proteins[END_REF], where the authors control the accessibility of the binding site with a square wave of unit amplitude h(t) that enters the boundary condition of the diffusion equation 4πR

Examples of time-dependent gating functions

Our goal is to develop an analysis of the gated reactions which involves the use of arbitrary shaped deterministic gating functions. This choice can be useful for modeling reactivity fluctuations in which the opening/closing happen gradually, and where the time of opening and closing (both the permanence in a given state and the velocity with which the accessibility changes) can be controlled by appropriate parameters.

Before introducing the two classes of function that we adopted to model the gating function, we want to stress that in order to solve the Volterra equation with the method described in the previous section we have to require the function h(t) to be continuous. Bearing in mind this restriction, we consider two kinds of periodic functions for the gating h(t): The rate constant of a reflecting sphere featuring an absorbing patch of aperture θ 0 has been computed in Ref. [START_REF] Traytak | Diffusion-controlled reaction rate to an active site[END_REF]. The result is k = f (θ 0 )k S (6.24) where the steric factor that accounts for the rate reduction can be well approximated for small angles θ 0 by: f (θ 0 ) ≈ θ 0 + sin θ 0 2π -(θ 0 + sin θ 0 ) . (6.25)

The analysis reported in Ref. [START_REF] Traytak | Diffusion-controlled reaction rate to an active site[END_REF] can be easily extended to the case of a partially reactive patch, that is, a patch connoted by an intrinsic reaction rate constant k * . In this case, one finds: where θ 0 and α are fit parameters.

k k S = H 0 f (θ 0 ) H 0 + f (θ 0 ) . ( 6
We conclude that the mean rate constant for a fully reactive spherical receptor featuring a gating dynamics can be mapped onto the rate constant of an equivalent static receptor featuring 1) a reduced reactive surface (patch) and 2) an effective (reduced) intrinsic rate constant of the reactive patch. Chapter 7

Conclusions and perspectives

The general purpose of this thesis is to analyze diffusion-reaction processes in different types of non-ideal conditions. Among them, we take into account the interaction among species that are diffusing in the same space, the effects of crowding on the motility of agents within densely populated media and the presence of obstacles, compartmentalization and traps in the domain. These are fundamental aspects for describing many biochemical processes in living media, as well as for studying industrial and chemical applications of diffusion influenced reactions. As underlined by several experimental results, when the features of the domain where the processes occur are highly non-ideal (high density, complex geometries) far from the ideal approximation, the classic diffusion theory which applies to diluted media is no longer adequate. Neglecting the interaction among diffusing agents and between the agents and the external environment prevents one to grasp the essential features of the processes of interest. This work aims at introducing some degree of complexity in the description of diffusion and diffusion-reaction processes, and at identifying to which extent these non-ideal conditions influence the mobility of particles and the rate of the reactions occurring between molecules.

This thesis can be divided into two main blocks. The first part is devoted to macromolecular crowding, namely the presence of a finite density of reactants and other inert species, which is a common feature to many realistic situations in biology and soft matter. For all the processes described in the first chapter, the crowding is enforced in the microscopic description of a random walk, by accounting for excluded volume effects in the analysis of the single agents' movement. This approach allows us in primis to recover the macroscopic description of the evolution of the particle density for different choices of the jump rules in the mean-field limit. We show that considerable information is lost in passing to the continuum: the excluded volume effect introduced at the microscopic level disappears unless the description accounts for some degree of asymmetry, such as in the presence of anisotropic domains or biased rules for the motion of the agents. Deviations from the standard diffusion equation emerge also by making part of the diffusing population distinguishable, or by taking into account explicitly the size of the agents in the extended crowding scenario. The microscopic models themselves can provide interesting information supplementary to the mean-field description, such as the role of the fluctuations. Such analysis is not among the purposes of this thesis, but could be carried out starting from the agent-based models.

In the second part of this work, the Smoluchowski theory for diffusion-influenced reactions is adapted to non-ideal domains, where many reactive boundaries interact with the diffusing ligands and compete for the same resources. The role of the environment's geometry (obstacles, compartmentalization) and distributed reactivity (competitive reactants, traps) is quantified through the explicit calculation of the reaction rate constant in a domain arbitrarily decorated with obstacles and reactive boundaries. The only constraint imposed is that all bodies have to be represented as assemblies of spheres. For example, complex-shaped molecules, such as antibodies, can be represented through coarse grained models that employ spherical beads as elementary components. The diffusion-reaction process, initially defined as a boundary value problem for the density of the ligand molecules, with boundary conditions enforced on the surface of the spheres, is mapped onto a linear system of algebraic equations through the use of harmonic expansion techniques. This procedure enables us to compute analytically the reaction rate constant of the selected reaction to any desired accuracy, with the employ of analytical methods only. We show that this approach provides remarkable advantages with respect to the numerical methods based on finite-element computations, because it leads to a much more accurate solution, and because it allows one to derive simple approximate analytical formulas. We show that this theoretical framework can be used to investigate naturally occurring binding processes and to assist in the design of artificial structure such as composite nano-reactors used in industrial catalysis. Moreover we demonstrate that in many cases, such as in the description of the antibody-antigen binding process, a complex system can be mapped onto a simpler schematic model that, even though with less detail, is still able to explain the basics of the process. We stress that our goal is to devise a mathematical framework to characterize the interplay of diffusion and reaction with static environmental factors. Obviously, even though we relaxe some of the basic hypotheses of the classical ideal description, we still work under simplified assumptions with respect to real media. Nevertheless, our approach allows us to identify important trends and correlations between physical-chemical and geometrical parameters, such as the configuration of the particles in a multi-catalyst nano-reactor that optimize the reaction. These aspects could not be recognized in the infinite-dilution treatment, which neglects the multi-body interactions. The procedure followed to solve the stationary diffusion equation in the presence of multiple boundaries could be in principle extended to different geometries realized with more complex, non-spherical bodies. The only requirement is the existence of addition theorems for the elementary solutions in the selected coordinate system, which would allow one to express the boundary conditions over the multi-body landscape.

As an example, we can imagine to apply the same method to geometries involving elongated agents, described in ellipsoidal coordinates.

Our method can be extended to the Laplace space, so as to work out exactly the effect of the environment on time-dependent problems. In this case, one has to employ addition theorems for the Bessel functions, and the difficulty lies in the subsequent inversion of the Laplace transform [START_REF] Sd Traytak | On the time-dependent diffusive interaction between stationary sinks[END_REF]. In our treatment of diffusioninduced reactions, the time-dependent effects are considered only with respect to a simplified geometry. This is the topic of the final chapter of this work, where we take into account a temporal variation of the reactivity of a single isolated receptor. The analysis is performed with a different method from the previous treatment, but with the same result of transforming a boundary-value problem into a simpler mathematical structure. The differential problem with complex time-dependent boundary conditions is in fact mapped onto an integral Volterra equations which can be solved through well-known and stable iterative procedures.

The method derived to characterize arbitrary shaped geometries can be used to investigate some geometric effects on chemoreception. Cellular receptors are not uniformly distributed. For example clusters of receptors are recognizable on selected regions of the surface of bacteria as in E.Coli, and their spatial arrangement resembles lattice-like structures. One can investigate the reason why some arrangements seem to be preferred to others, and relate them to the optimized configuration. Along these lines, we are currently studying the effect of number, dimension and configuration of the receptors on the overall reaction rate to a cell surface. The receptors can be represented as spherical absorbing boundaries which cover part of the surface of an otherwise reflecting sphere (see Fig. 7.1). Our preliminary results indicate that ordered structures of receptors maximize the overall rate to capture, while extended localized receptor-covered patches correspond to very low capture rates. These preliminary results point to a highly non-trivial connection between the geometry and the physics of chemoreception.

Marta GALANTI

Processus de diffusion et réaction dans des milieux complexes et encombrés L'objectif général de cette thèse est d'analyser les processus de diffusion et les processus de réactiondiffusion dans plusieurs types de conditions non-idéales, et d'identifier dans quelle mesure ces conditions non idéales influencent la mobilité des particules et les réactions entre les molécules. Dans la première partie de la thèse, nous nous concentrons sur les effets de l'encombrement macromoléculaire sur la mobilité, ainsi élaborant une description des processus de diffusion dans des milieux densément peuplés. Tous les processus sont analysés à partir de la description microscopique du mouvement des agents individuels sous forme de marche aléatoire, tenant compte de l'espace occupé par les particules voisines. La deuxième partie de la thèse vise à caractériser le rôle de la géométrie de l'environnement et de la réactivité des corps qui y sont contenus sur la réaction entre des molécules sélectionnées. La théorie classique de Smoluchowski, formulée pour les réactions contrôlées par la diffusion dans un milieu dilué, est ainsi adaptée à des domaines arbitrairement décorés par des obstacles, dont certains réactifs, et l'équation stationnaire de diffusion est résolue avec des techn Finalement, le calcul explicit de la constante de réaction et la dérivation des formules approximées sont utilisés pour étudier des applications biologiques et nano-technologiques.

Mots-clés : diffusion, réaction-diffusion, crowding moléculaire, Smoluchowski, analyse harmonique.

Diffusion-reaction processes in complex and crowded environments

The overall purpose of this thesis is to analyze diffusion processes and diffusion-reaction processes in different types of non-ideal conditions, and to identify to which extent these non-ideal conditions influence the mobility of particles and the rate of the reactions occurring between molecules. In the first part of the thesis we concentrate on the effects of macromolecular crowding on the mobility of the agents, providing therefore a description of various diffusion processes in densely populated media. All the processes are analyzed by modeling the dynamics of the single agents as microscopic stochastic processes that keep track of the macromolecular crowding. The second part of the thesis aims at characterizing the role of the competitive reactants, traps) on the reaction between selected molecules. The Smoluchowski theory for diffusioninfluenced reactions is thus adapted to domains arbitrarily decorated with obstacles and reactive boundaries, and the stationary diffusion equation is explicitly solved through harmonic-based techniques. The explicit calculation of the reaction rate constant and the derivation of simple approximated formulas are used for investigating nano-technological applications and naturally occurring reactions.

Figure 2 . 2 :

 22 Figure 2.2: Subsequent frames in a typical Single Particle Tracking experiment.

Figure 2 . 3 :

 23 Figure 2.3: Fluorescence intensity of the monitored volume against time. The steady state intensity differs from the initial value because of the fraction of molecules which are immobilized at the bleached spot.

Figure

  Figure 2.4: Comparison between stochastic (black solid lines) and mean field simulations (red solid line). The stochastic simulations are averaged over 20000 realizations. Here φ(x, 0) = 0.5. Two dimensional snapshots of the stochastic simulations are also displayed. The crowders are plotted as small black circles. The tagged particle is represented by the filled red square. Upper panel: initial condition (see first snapshot, top-left panel) originating a super-diffusive transient. Lower panel: initial condition (see first snapshot, top-left panel) causing a sub-diffusive transient. The two kinds of initial condition are the same considered in section 2.2.6

Figure 2 . 5 :

 25 Figure 2.5: Rescaled MSD of the tagged species as function of time for different crowding strengths in 1D. Time is expressed in units of τ 1 = x 2 0 /2D.

Figure 2 . 6 :

 26 Figure 2.6: Rescaled MSD of the tagged species as function of time for different crowding strengths in 2D. Time is expressed in units of τ 2 = r 2 0 /2D. The qualitative behavior is the same observed for the 1D domain.

Figure 2 . 7 :

 27 Figure 2.7: Mean square displacement for a pulse of the kind (2.49) with β = 10 as a function of time for different choices of the crowding parameter φ. The dashed lines are straight lines of slope v e as predicted by eq.(2.48). Inset: approach to the diffusive regime. The parameter D = v 2 /2r is the theoretical diffusion coefficient. The solid line refers to the PRW and vanishes as t -1 . The dashed line is an inversepower law with exponent 0.4, to be used as a guide for the eye.

Figure 2 . 8 :

 28 Figure 2.8: Effective velocity (2.50) as a function of the intensity of crowding for a Gaussian pulse , for a generalized Gaussian with β = 10 and for a step pulse (β → ∞).

Figure 3 . 1 :

 31 Figure 3.1: Scheme of the interior of a cell showing the high density and heterogeneity of cytoskeletal structures.

Figure 3 . 2 :

 32 Figure 3.2: Scheme of the bimolecular reaction (3.1).

Figure 3 . 3 :

 33 Figure 3.3: The problem of Smoluchowski

Figure 3 . 4 :

 34 Figure 3.4: Schematic representation of the domain Ω = Ω 0 \

Figure 3 . 5 :

 35 Figure 3.5: Multiple coordinate system for a general domain.

Figure 3 . 6 :

 36 Figure 3.6: Scheme of the the application of the addition theorems (3.10), (3.22) and (3.24) for expressing the boundary conditions in the local coordinates on ∂Ω 0 and on each ∂Ω α . The choice of the appropriate addition theorem depends on the ratio between the distance between the centers of each pair of reference systems and the norm of the position vector in the new reference system (see the discussion in subsection 3.1.1.)

Figure 3 . 7 :

 37 Figure 3.7: Vectors expressing the position of a point P in two different coordinate system S α and S β

. 38 )

 38 Conversely, if |r α | > |L βα |, then we have to write the irregular harmonic as a function of the irregular harmonics centered in S α :

Figure 3 . 8 :(- 1 )

 381 Figure 3.8: Schematic representation of the domain R 3 \ Ω 1 where we look for the solution of the boundary problem (3.40).

. 47 )Figure 3 . 9 :

 4739 Figure 3.9: Schematic illustration of an axially symmetric arrangement of spherical boundaries Ω α inside a larger spherical domain Ω 0 .

Figure 3 . 10 :

 310 Figure 3.10: Rate constant normalized to the rate of an isolated sink of the same size for a two-sink configuration with R α /R β = 5 in the MOA (3.55). The rates are plotted versus the distance L/R

Figure 3 . 11 :

 311 Figure 3.11: Normalized rate to the large sink in a system of two spheres for different choices of the ratio c = R 1 /R 2 . The rate is plotted versus the distance δ between the boundaries, normalized to the small sink's linear dimension. The symbols are the results of the analytical computations, the solid lines are plots of eq. (3.55). The position of the minimum, that is absent for c = 1, is shifted to the right as c is increased.

Figure 3 . 12 :

 312 Figure 3.12: Schematic illustration of a sink Ω 1 confined within a spherical cavity Ω 0 with the centers of the two spheres lying on the z-axis at a distance L.

Figure 3 . 13 :

 313 Figure 3.13: Two reactive spherical boundaries inside a spherical cavity separated from the outside by a permeable barrier.

Figure 3 . 14 :

 314 Figure 3.14: Diffusion interaction strength against the normalized distance between two identical sinks of radius R. The computation (symbols) is performed for two values of R. The monopole (dotted line) and dipole approximation (solid line) are compared to the exact results.

Figure 3 . 15 :

 315 Figure 3.15: Diffusion interaction strength against the position of the barycenter of the system of two sinks, for two different choices of the distance L. The variables z 1 and z 2 are the z-coordinates of the centers of Ω 1 and Ω 2 .

Figure 3 . 16 :

 316 Figure 3.16: Total ligand flux to the central sink in a two-neighbor model normalized to k S = 4πDσ. Two cases are shown for the neighbors, reactive (red) and reflecting (blu boundaries).The rate is plotted as a function of the normalized distance to the central particle. The analytical solution of (3.67) is compared to the numerical results taken from Ref.[START_REF] Eun | Influence of neighboring reactive particles on diffusion-limited reactions[END_REF]. The different symbols correspond to different angles between the neighbors.

  ο , FE num. sol.

Figure 3 . 17 :

 317 Figure 3.17: Comparison between the numerical solution [117], our solution, and different degrees of approximation in the perturbative expansion (up to the first, second and third order in ε) for N = 2 neighboring spheres separated by an angle of 60 • . One can remark the alternation in sign of the consecutive terms of the expansion. The rate constant is represented versus the distance to the central sink, normalized to the diameter of the particles.

Figure 3 . 18 :

 318 Figure 3.18: Perturbative (eq. (3.71)) vs exact rates to a sink of radius σ screened by N = 50 sinks of radius σ 1 = σ/10 arranged randomly at a fixed distance d from it. Each set comprises 100 independent configurations. The values of the rate are normalized to the Smoluchowski rate for the central sink, k S = 4πDσ. the perturbative rate constants are to an excellent extent proportional to the exact values.

. 72 )

 72 with P (ω αβ ) = sin ω αβ /2, and 2 arcsin(σ 1 /d) ≤ ω αβ ≤ π (excluded-volume constraint between screening sinks). The probability density is normalized such that

Figure 3 . 19 :Fig. 3 .

 3193 Figure 3.19: Total flux into a sink of radius σ surrounded by N spherical boundaries of the same size arranged randomly at distance d (normalized to k S = 4πDσ). Symbols denote the exact results (solution of the linear system), each point marking the average rate computed over 100 independent configurations, while the shaded bands highlight the regions comprised between the minimum and maximum rates. For reflecting screening boundaries, these regions are as small as the truncation error. The light blue and orange lines are plots of the perturbative expansion for the inert neighbors(3.70) and of the monopole approximation (3.73) for the reactive boundaries, respectively, with λ = 1. Two configurations with N = 50 are shown explicitly with the screening spheres made all absorbing (bottom) and all reflecting (top), with arrows flagging the corresponding values of d and normalized rates. The figure illustrates that the competitive screening greatly reduces the rate constant compared to inert obstacles, and it is strongly modulated by the configuration.

75 )

 75 Our goal is to obtain an expansion of the rate to the central sink, k/k S = -B 0 00 in powers of the parameter ε = σ/d. For this purpose, we have to consider the

Figure 3 . 20 :

 320 Figure 3.20: Normalized rate to a central sink surrounded by N = 34 spheres, for inert and reactive boundaries. Each symbol represents the result corresponding to a given configuration, and the arrangements are ordered on the base of the standard deviation of the distances between the spheres in a given conformation.

Figure 3 .

 3 Figure 3.21: Configuration-dependent parameters in the perturbative expansion (3.71).

Figure 3 . 22 :

 322 Figure 3.22: Total flux into a sink of radius σ surrounded by N = 50 smaller sinks of radius σ 1 = σ/10 arranged randomly at a distance d (normalized to k S = 4πDσ). The left-most and right-most cartoons depict two configurations that screen exactly the same amount of flux, despite being at considerably different distances (d/σ = 1.1 and d/σ = 8). The configuration shown in the middle corresponds to the predicted minimum at d/σ = 1 + 1 -λ[1 + (N -1)λ] ≈ 1.64. The solid line is a plot of formula (3.73). Each symbol is the average over 250 independent configurations, while the filled band comprises the region between the minimum and maximum rates. The top panel illustrates the case of screening by a large number of tiny particles, highlighting the sizeable non-monotonic effect. The curves are plots of eq. (3.73).

Figure 4 . 1 :

 41 Figure 4.1: Scheme of a core-shell nanoreactor with embedded gold nanoparticles. The PS core is shown at the center.

Figure 4 . 2 :

 42 Figure 4.2: Schematic representaation of a composite nanoreactor displaying the parameters used in the analytical derivation.

  denotes the interior of the PS core and Ω α = {|r α | = |r 0 -L α | < a}, α = 2, 3, . . . , N + 1, denote the interior of the α-th nanosphere. The ligand diffuses with diffusion coefficients D i and D o inside the microgel shell and in the bulk, respectively (D i < D o ).

Figure 4 . 3 :

 43 Figure 4.3: Reaction rate constant versus number of nanoparticles of radius a = 0.0146R 0 for D i = D o (upper panel, swallen phase) and D i = 0.2D o (lower panel, shrunk phase). Symbols denote the exact results, while the solid lines are plots of the configuration-averaged monopole approximation, eq. (4.28). Other parameters are R S /R 0 = 0.353, h → ∞.

Figure 5 . 1 :

 51 Figure 5.1: Schematic representation of the polypeptidic chains forming an antibody molecule. The yellow bonds show explicitly the structure-stabilizing disulfide bridges.

Figure 5 . 2 :

 52 Figure 5.2: Pictorial representation of the three-lobe structure of an antibody. VL and VH denote the highly variable regions of the Fab domains where antigens bind.

Figure 5 . 3 :

 53 Figure 5.3: Coarse grained representation of the IgG with N = 96 beads. The docking of the N spheres is performed with the shape-based algorithm[START_REF] Arkhipov | Stability and Dynamics of Virus Capsids Described by Coarse-Grained Modeling[END_REF][START_REF] Arkhipov | Four-scale description of membrane sculpting by BAR domains[END_REF] implemented in the VMD package[START_REF] Humphrey | VMD -Visual Molecular Dynamics[END_REF]. The black spheres positioned on the edges of the Fab represent the active sites. In this schematic illustration, the beads occupy the positions that are assigned through docking. The radius of the beads shown in the figure corresponds approximately to the one used in the simulations (0.44 nm). The atomic representation of the antibody is shown as a transparent structure (PDB code: 1IGT)

Figure 5 . 4 :

 54 Figure 5.4: Snapshots obtained from the coarse-grained molecular dynamics simulations. The points represent the position of the centers of the units of the coarsegraining. The dimension of the spheres is purely indicative of the position and is not representative of the radius used in our model.

Figure 5 . 5 :

 55 Figure 5.5: (a) Cryo-ET procedure: single molecules in vivo are immobilized at cryogenic temperatures (below -150 • C) and independent snapshots of different sections of the body are taken by rotating the sample inside the cryo-microscope. The images are then processed and assembled to reconstruct the three-dimensional structure of the molecule. (b) A set of reconstructed 3D volumes of the IgG molecules in solution. The three different tomograms representing individual IgG molecules.The image is reproduced from Ref.[START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF].

0 ∂u ∂r α ∂Ωα = 0 for α = 1 , 2 u|Figure 5 . 6 :

 1256 Figure 5.6: (a) Cumulative distribution of the Fab-Fc angle φ. The shaded area identify the regions that are inaccessible because of steric hindrance. (b) Cumulative distribution of the Fab-Fab angle ψ.The solid line stands for the results of our coarse-grained simulations, while symbols refer to the Cryo-ET experiments reported in Ref.[START_REF] Bongini | Freezing immunoglobulins to see them move[END_REF]. The inset in panel (a) provides a schematic view of the mechanical model employed in the simulation, together with a definition of the relevant angles ψ and φ.

Figure 5 . 7 :

 57 Figure 5.7: Schematic illustration of the set of boundaries corresponding to a given conformation of an antibody molecule. Blue spheres correspond to reflecting surfaces, while red beads flag absorbing boundaries (paratopes).

Figure 5 . 8 :

 58 Figure 5.8: Normalized rate against the distance between the active sites (red beads in the snapshots) for two different choices of the encounter distance R a = paratope size + antigen size. The rate constant is normalized to twice the Smoluchowski rate constant of an isolated paratope, k S = 4πDR a , where D = D IgG + D A ≃ D A is the relative diffusion constant (practically equal to the antigen diffusion coefficient). Each symbol represents the rate calculated for a given configuration of the IgG. The radius of the reflecting spheres that define the body of the IgG was fixed at 0.44 nm in both cases. The solid and dotted lines refer to the effective model 3.53 multiplied by a constant factor f a . Here f a = 0.94 (R a = 0.54 nm) and f a = 0.93 (R a = 0.44 nm). All the calculations are performed with N = 96.

Figure 5 . 9 :

 59 Figure 5.9: Histograms of the reaction rate constants computed for different choices of the paratopes' size R a . The solid lines are fits to the Gumbel density distribution P(x) = exp (xµ)/βe (x-µ)/β /β, where µ, β are the fitting parameters. The best-fit values of the parameters are: µ = 1.81, β = 0.010 (R a = 0.44 nm), µ = 1.77, β = 0.012 (R a = 0.54 nm).

  Figure 5.10: Left: rate constant normalized by 2k S versus rescaled paratopeparatope distance for different choices of the paratope size R a . The grey dashed lines are plots of formula (3.53), rescaled by appropriate normalization factors f a as shown in the right panel (the solid line is only a guide to the eye). The black dotted line is the plot of the perfect dumbbell prediction (f a = 1), k/k S = d/(d + R a ).

  on arbitrary placed boundaries ∂Ω i in the open domain Ω = R 3 \ N i=1 Ω i . More specifically Duhamel's theorem provides a solution for the problem

Figure 6 . 1 :

 61 Figure 6.1: Average rate against the maximum amplitude of the gating function h(t) = H 0 (cos(ωt) + 1)/2 for different choices of ω. The initial value is ρ 0 (r) = (1 -1/r). .

. 26 )

 26 In fig6.3 we illustrate the fits of the average rate constant corresponding to symmetric and asymmetric gating functions to expression (6.26) with an effective reactivity αH 0 :k k S = αH 0 f (θ 0 ) αH 0 + f (θ 0 ) . (6.27)

Figure 6 . 2 :

 62 Figure 6.2: (a)Plot of the square wave gating function with asymmetric opening and closing time constants.(b) Mean rate for different choices of the closing time ∆t 2 . Other parameters are τ O = 8τ D , τ C = 5τ D , ∆t 1 = 0.1τ D and the integration time is T = 70τ D The initial value is ρ 0 (r) = (1 -1/r).

Figure 6 . 3 :

 63 Figure 6.3: Average reaction rate versus the maximum amplitude of the gating wave function for symmetric and asymmetric nearly square-wave gating functions of the form(6.22). The data are fitted to expression (6.27) via the parameters θ 0 and α. The choices of the parameters in the gating functions are τ O = 8, τ C = 5, ∆t 1 = 0.1. The integration time is T = 70τ D The initial value is ρ 0 (r) = (1 -1/r).

Figure 7 . 1 :

 71 Figure 7.1: Analysis of the total ligand flux to N = 100 receptors of radius a = 0.03 R on the surface of a reflecting sphere of radius R. The histogram shows the distribution of the rate over a population of over 2 × 10 4 independent random configurations of the receptors (normalized to the population average). Explicitly shown are a) the configuration corresponding to the least possible flux (all receptors in a spherical patch), b,c) the configurations corresponding to the least and the largest rates over the random ensemble, d) the configuration corresponding to the largest possible rate. Interestingly, the latter is identical to the configuration that minimizes the Coulomb energy of N identical charges on the sphere (Thomson configuration).

  Key words: diffusion, reaction-diffusion, crowding, Smoluchowski theory, harmonic analysis Università degli Studi di Firenze Via di Santa Marta 3, 50139 Firenze

  .66) Therefore, recalling eqs. (2.61) and (2.62), we have from eqs. (2.63), (2.64) and (2.65), respectively,

Table 2 .

 2 1: The coarse grained particle conditional pair distribution function Ḡp for stationary ensembles of d-dimensional spheres of diameter σ[START_REF] Torquato | Nearest-neighbor distribution functions in many-body systems[END_REF]

  2 D ∂ρ(r, t) ∂r | R = h(t)k * ρ(R, t),(6.17) which has the same structure as the BC in our system (6.1). The equation is solved by using finite difference numerical methods for solving the equivalent formulation reflecting boundary condition at r = R and with an outer boundary condition ρ(∞, t) = 1. Once the boundary value problem is solved numerically for the density of ligands, the authors compute the time dependent rate coefficient k(t) corresponding to a gate which is alternatively open and closed for intervals of arbitrary length τ O and τ C . They provide an analytical approximation for k(t) in terms of the ungated time-dependent rate k u (t) which holds in the limit τ C ≫ τ D , where τ C and τ D stand for the gate closing time and the relaxation time, respectively. This assumption allows the density of ligand to recover the initial value ρ(x, 0) during the closing interval. The average rate constant is evaluated over an open/closed cycle as the average rate constant of the ungated case times the fraction of time that the gate is open: Two limiting behaviors within the regime τ C ≫ τ D are retrieved for short and long opening times. If τ O ≫ τ D the time-averaged gated constant is equal to the fraction of time the gate is open, h(t) , times the stationary limit of the rate for the ungated binding:

	∂ρ ∂t	= r -2 ∂ ∂r	r 2 D	∂ρ ∂r	-	k * 4πR 2 h(t)ρδ(r -R)	(6.18)
	coupled with a k = h(t)	1 τ O	0	τ O	k u (t)dt =	1 τ O + τ C	0	τ O	k u (t)dt.	(6.19)
				k ∼	k * k D k * + k D	h(t) .	(6.20)
	If τ

O ≪ τ D the average normalized gated constant becomes k ∼ k * h(t) .

(6.21)

We emphasize that we use the general terminology of d-dimensional hard spheres. Obviously, these are hard rods in one dimension and hard disks in two.

We note that the same symbol kS is used in the literature to refer to the rate constant, with dimension [M -1 s -1 ]

(a) (b)

Conclusions and perspectives

Bibliography

for n = 0 and n = 1, namely

.

(3.77) Substituting the expressions (3.76) in the above equations, the unknown coefficients Q α,j m,n are determined by equating the coefficients of equal powers of ε. With the use of basic trigonometric identities, including θ 0β = πθ β0 and φ 0β = π + φ β0 , we get the first coefficients of the expansion:

• a continuous periodic function of arbitrary amplitude H 0 that connects the open and closed states, represented by intervals of arbitrary length τ O and τ C connected by raising stages of duration ∆t

By prescribing very steep rises, this function mimics a square wave (the function used in Ref. [START_REF] S H Northrup | Rate theory for gated diffusioninfluenced ligand binding to proteins[END_REF]). The parameter H 0 represents the reactivity of the receptor when the gate is completely open.

• A periodic wave of amplitude H 0

The frequency ω is the parameter which controls the fast/slow gating regimes.

For all the choices of the functions and parameters we computed the time dependent rate constant as well as its mean value in the integration time. In general, for a periodic gating function, the integration time T is selected as a multiple of the period so that the running average reaches a stable value. In our analysis we considered D = 1, R = 1 and ρ B = 1.

In the slow gating limit we found the same value for the averaged rate constant obtained in the previous studies. For the fast gating regime, we could not compare directly our findings with the limit recovered in the stochastic analysis, because the use of a continuous gate function does not allow to model a kinetics where the fluctuations happen infinitely fast. Nevertheless we model this particular regime with a periodic sine wave (the second class of functions introduced before) with high frequency. In figure 6.1 we represent the average rate constant as a function of the amplitude of h(t) for different choices of the frequency. We can see that the qualitative behavior of the rate as a function of the maximum reactivity is the same as the ungated case, and the average rate tends to the corresponding stationary ungated limit as ω is increased.

An interesting phenomenon that we want to investigate is the correlation between symmetry in the characteristic open/closing time and binding efficiency. More specifically we wonder how the rate constant is modified by acting on the parameters which govern the opening and closing of the gate. In figure 6.2 we plot the average rate as a function of the maximum amplitude of a modified gating function of the kind (6.22), where the closing and opening dynamics are described by two different time constants ∆t 1 and ∆t 2 (see Fig 6 .2, top panel), and we compare the behavior for different ratios of the rising stage parameters ∆t 1 and ∆t 2 . One notices immediately that the binding efficiency is increased as the symmetry of the function which regulates the gating of the site is lost. We stress that the average reactivity (which is given by the integral of h(t) over an arbitrary number of periods, divided by the time of integration) is the same for all the considered gating functions h(t). We can thus appreciate here how the asymmetry enhances functionality.

The kinetic analysis of the gated reaction rate to a homogeneously reactive sphere can be mapped onto the stationary binding with a spatially anisotropic reactant.