Simulation-based Reconstructed Diffusion unveils the effect of aging on protein diffusion in Escherichia coli

We have developed Simulation-based Reconstructed Diffusion (SbRD) to determine diffusion coefficients corrected for confinement effects and for the bias introduced by two-dimensional models describing a three-dimensional motion. We validate the method on simulated diffusion data in three-dimensional cell-shaped compartments. We use SbRD, combined with a new cell detection method, to determine the diffusion coefficients of a set of native proteins in Escherichia coli. We observe slower diffusion at the cell poles than in the nucleoid region of exponentially growing cells, which is independent of the presence of polysomes. Furthermore, we show that the newly formed pole of dividing cells exhibits a faster diffusion than the old one. We hypothesize that the observed slowdown at the cell poles is caused by the accumulation of aggregated or damaged proteins, and that the effect is asymmetric due to cell aging.


Introduction
Diffusion of molecules inside cells plays a crucial role in the func1oning of biochemical processes. In prokaryo'c cells, which lack membrane-bound compartments, except for the periplasm in Gram-nega%ve bacteria, the majority of cellular processes take place in the cytoplasm. Here, the random mo#on of (macro)molecules allows for the func#oning of highly sophis#cated systems, such as the correct localiza+on of the septa+on ring in Escherichia coli, governed by a reac.on-diffusion mechanism 1 , the signal transduc"on that leads to chemotaxis 2 , or the transla,on of mRNA into proteins by polysomes, ribosomes-mRNA assemblies, mostly localized at the cell poles [3][4][5][6] .
The cytoplasm of bacteria is an extremely crowded environment, where the concentra7on of macromolecules, mainly proteins and RNAs, can reach values up to a volume frac8on of 15-20% in growing cells [7][8][9] and even higher in osmo(cally stressed cells [10][11][12][13] . The molecular composi0on of the cytoplasm is highly diverse, with molecules spanning in size over more than three orders of magnitudes, from sub-nanometric for ions and metabolites, to micrometric for the chromosome 14,15 . Despite this varia%on in size and surface proper%es of the molecules, many cellular components are uniformly distributed throughout the cell 5,16,17 , with notable excep0ons like the chromosome and nucleoid-binding proteins (localized in the cell center 4,18 ), the polysomes (localized at the poles and cytoplasmic periphery 3,4,6 ), and aggregated or misfolded proteins (localized at the cell poles [19][20][21] (Fig. 1). In addi-on, there is increasing evidence for the forma1on of phase-separated liquid droplets or biomolecular condensates in the cytoplasm of microorganisms [22][23][24] , which are metastable structures where certain proteins par33on.
Diffusion of spherical par1cles in aqueous solu1ons can be described by the Einstein-Stokes equa#on 25 . However, mo*on of par*cles in the highly crowded and inhomogeneous cytoplasm of bacterial cells has been extensively documented to deviate from the aforemen5oned model 16,17,26,27 . We have shown that the apparent diffusion is solely dependent on the complex mass, that is the molecular weight of the monomer, summed with the molecular weight of the fluorescent reporter, mul6plied by the oligomeric state of the analyzed molecules 16 (Fig. 1). The same conclusion was recently obtained in another study, using a different method 17 . Apparent diffusion of proteins interac2ng with large cellular components, such as the cell membrane, the nucleoid, the ribosome or the proteasome complex, depends on the sum of both the interac-ng masses and on the interac-on strength (dissocia-on constant, K D ) of the species 17,28 . The devia*on from the Einstein-Stokes equa+on indicates that diffusion in cells depends not only on the size of the analyzed molecules, but also on the composi5on and physical state of the cytoplasm. This, in turn, is dependent of e.g. (fluidiza4on by) metabolism 29,30 , catalysis-induced enzyme movement 31,32 and environmental stresses 11,33 . The devia*on from Stokes-Einstein can be explained by taking into account both the solvent quality of the cytoplasm and the size of the diffusing species. The macromolecular crowding experienced by each protein is dependent of its own molecular weight, that is, smaller molecules will be less affected by the crowded cytoplasmic environment than bigger ones 16,17,26,27,34 (Fig. 1). We have shown that this so-called perceived macromolecular viscosity is not spa2ally uniform in the cell 16 .
Despite the advancements made in the field of single molecule fluorescence microscopy 16,17,[35][36][37][38][39] , diffusion measurements are s.ll highly influenced by the effect of confinement, especially in small compartments such as the bacterial cytoplasm 16 and periplasm 40,41 , and eukaryo,c organelles 42 . Here, diffusion coefficients near the cell boundaries always appear lower than in the cell center 16,38 . This does not allow to properly separate the effect of confinement from possible physiological slowdown in diffusion of the analyzed species 16 (Fig. 1). Moreover, techniques such as Single Par:cle Tracking (SPT) and Single Molecule displacement Mapping (SMdM) produce a two-dimensional output of a threedimensional mo+on, which leads to obvious shortcomings in the es+ma+on of diffusion coefficients of par$cles freely moving in the cell cytoplasm.
Some methods to resolve confined diffusion have been developed. Bickel 43 proposed a mathema,cal method to obtain the mean square displacement in disks and sphere for par5cles. Bello8o et al. 17 derived a Ornstein-Uhlenbeck model for fi0ng Fluorescence Correla7on Spectroscopy (FCS) data acquired in a confined cylinder of infinite length, but not considering the effect of the cell poles, which represent the zones where the confinement affects the observed diffusion the most 16 .
A detailed analysis of the diffusion of macromolecules in the cell poles is paramount to understand the effects of aging in rod shaped cells. A study that followed repeated divisions of E. coli suggests that cells that inherit the old pole exhibit a diminished growth rate, decreased offspring produc7on, and an increased incidence of death 44 . Asymmetry in the doubling 3me of old and new pole daughter cells has also been observed in a more recent study in E. coli 45 . The underlying mechanisms of aging of bacteria are at best poorly understood.

Fig. 1. Diffusion in the cytoplasm of E. coli.
A diffusion map obtained with SMdM is overlayed with a schema9c of the cytoplasm of the cell. The top panel highlights the effect of confinement on the measured diffusion, which leads to lower diffusion coefficients near the boundaries of the cell. The bo"om panel represents the effect of the perceived viscosity by diffusing proteins. Since diffusion scales with the complex mass, bigger par$cles will be affected more by the crowding of the cytoplasm than smaller molecules and move rela$vely more slowly, leading to the devia.on from the Einstein-Stokes equa+on. The le1 panel represents our current hypothesis on the observed slowdown at the cell poles compared to the cell center, with accumula6on of aggregated or misfolded proteins impairing the diffusion in these regions.
In this study, we set out to solve the shortcomings of previous methods for analyzing confined diffusion in small compartments, and we apply the new tools to inves4gate the rela4onship between diffusion in the cell pole regions and aging. We developed a method to determine diffusion coefficients that are not affected by the effect of confinement or mo1on along the third dimension. We developed the technique to unveil novel proper,es of diffusion in small compartments and we demonstrate its validity in E. coli cells. We first inves.gate a mathema.cal approach, followed by a simula.on-based analysis, which we named Simula+on-based Reconstructed Diffusion (SbRD). We highlight the limita:on of using mathema&cal models to solve the problem of confinement, and the advantages of a simula7on-based approach. We use the method, in combina4on with a new cell detec4on tool, to re-analyze a previously acquired dataset 16 and obtain confinement-corrected values for the diffusion of molecules in the E. coli cytoplasm. Further, we inves6gate diffusion near the cell boundaries and at the cell poles to determine how much confinement influences the slowdown in these regions. We test the effect of an8bio8cs that disrupt the polysome structure in the apparent slowdown of diffusion at the cell poles. Finally, we make observa(ons about asymmetry in diffusion in the bacterial cell, which we associate with aging.

A mathematical solution to the limitations of confined diffusion
We previously showed through simula4ons that the measured diffusion coefficient is underes4mated when par)cles move by random mo)on in a confined environment. The devia)on from the real diffusion coefficient dras-cally increases when the diffusion coefficient gets higher, and the analyzed molecules are in a region closer to the boundary of the confinement 16 ( Fig. 2A).
We now developed a mathema/cal model to solve this limita/on. The simplest approach to modeling of diffusion is through Ficks laws, which in his second law led to the so-called diffusion equa.on (Eq. 1) In the context of random mo/on, say a random walk of a par/cle, ρ denotes a probability density (or distribu(on) and = d ! would be the probability of finding the par4cle within a region R.
Here we solve analy-cally the one-dimensional diffusion equa.on with reflec.ng boundaries. For simplicity, we interpret ρ(x, t) as the probability density for finding a par3cle at posi%on x at $me t.
Assuming that at t = 0 the par(cle is located at x = 0, and that the probability density ρ(x, t) = 0 as x approaches infinite for any finite 0me, a solu0on to equa0on 1 is then given by equa0on 2: In a heuris+c way, equa+on 2 provides a measure of the likelihood of finding a par+cle at posi+on x a"er !me t provided that the medium where it moves is infinite.
Let us now consider the diffusion of a par3cle within a bounded domain, for simplicity from -L/2 to + L/2, where L is the length of the domain. We assume that the par-cle is located at x = 0 at t = 0, and that it is reflected back to the interior of the interval once it reaches the boundary (Fig. 2B). The analy=cal solu%on for this case is (see Supplementary Informa5on -Diffusion on a closed interval): Relying on an analy+cal solu+on of the diffusion equa+on in higher dimensions and for complicated geometries is not convenient. Equa1on 2 is valid for diffusion in unbounded domains. We can approximate the solu/on of equa/on 3 using a "folding approach", by accoun/ng for the bounces a par$cle makes when hi0ng the boundaries, assuming no loss of energy in the process, and adding them up (Fig. 2C).
In this way, an approxima1on of equa1on 3 is given by equa1on 4: where the term ρ(kL -x) + ρ(kL + x) corresponds to the density for the par/cle being at x at $me t a"er k bounces. When subs,tu,ng equa,on 2 in equa,on 4 we get: where N denotes the maximum number of bounces (Fig. 2C).
One can take a similar approach in higher dimensions. For example, let us consider the diffusion equa%on in a rectangular plate of sides A (horizontal) and B (ver%cal), and let us assume that the mo%on of the par)cle on each coordinate (x, y) is independent of each other. Then, the analy4cal solu4on for the diffusion equa.on is: where ρ A and ρ B are as in equa*on 3.
We can however take another approach, as we did for the interval in one dimension (equa4on 4, 5). We denote by p(t) = (x(t), y(t)) the posi)on of a par)cle in the rectangular plate at )me t. As above, we can compute the density ρ(x, y, t) by adding up all the densi0es corresponding to trajectories in the rectangular plate that take the par.cle from the ini.al posi.on p 0 to some final posi-on p f a"er a number of bounces (Supplementary fig. 1A, 1B). What we describe is an example of a mathemaAcal billiard.
To use these ideas to es,mate the diffusion coefficient inside a cell, we approximate the geometry of the cell by a planar sphero-cylinder, also known as a Bunimovich stadium (Fig. 2D). In this se>ng, 0-bounce and 1-bounce trajectories can be easily computed. Densi5es corresponding to trajectories that bounce on the straight (top and bo0om) sides of the sphero-cylinder are computed as aforemen2oned. Those corresponding to bounces on the circular sides can be computed by solving the system of equa7ons: , (x f , y f ) represent the star*ng point, a bouncing point on the circular sec2on of the cell boundary and the end point, respec2vely (Fig.  2D).
Since equa*on 7 cannot be solved analy*cally, we developed a script (see Supplementary Materialalgorithm 1) to compute the 0-bounce and 1-bounce trajectories of any diffusing par4cle, provided that the start and end posi-ons are known. We then used this algorithm to calculate the diffusion coefficient of two-dimensional diffusion simula-ons in a billiard, generated with Smoldyn 46 (Fig. 2E). Smoldyn allows simula*on of the mo*on of par*cles using a predefined diffusion coefficient and *me resolu*on, within a simula,on compartment. With our mathema,cal model we obtain a final diffusion coefficient for every pixel that is a good approxima3on of the input diffusion coefficient used for the simula3ons. par$cle's posi$on at x(t). The first scenario is that the par$cle travels to its measured posi$on without bouncing. Another scenario is that the par2cle arrives at its measured posi2on a7er bouncing (on the boundaries) once. In this way, there are infinitely many ways in which the par2cle can reach its measured posi%on, depending on the number of bounces the par%cle made. However, the probability of each case is inversely propor,onal to the total distance traveled. In B we show the distances traveled for 0 and 1 bounce

Limitations of a mathematical approach
Applying the mathema.cal approach to solve the diffusion equa.on near the boundaries in confined environments has three shortcomings. (i) The approxima4on (equa4on 5) is valid for the given boundary condi&ons, but not for e.g. non-con$nuous, non-convex surfaces. An example of a surface where our model would have failed is the Penrose unilluminable room 47 (Fig. 3A), or the matrix of mitochondria. Some regions of these surfaces are inaccessible by rays that start from par5cular loca"ons, regardless of the number of bounces. However, for a par5cle freely diffusing in any compartment, it would be possible to reach any loca,on, leading to the emergence of star,ng and final points that cannot be connected by reflec%ng rays. (ii) The model cannot be extended from the two-dimensional billiard to the threedimensional spherocylinder. Given a start and end point in two dimensions, it is always possible to find a reflec%on point on a circle; on the other hand, given a start and an end point in three dimensions, there will be an infinite number of reflec2on points on a sphere. Therefore our model implies that the mo2on of par'cles only occurs in two dimensions (x,y coordinates of the diffusing par2cles), while in reality (in cells) par"cles also diffuse along the z-axis. When simula.ng diffusion in a three-dimensional spherocylinder and analyzing it with our mathema5cal model (equa5on 5), we observed an underes'ma'on of the diffusion coefficient throughout most of the cell, and an overes)ma)on of it close to the boundary (Fig. 3C). The underes9ma9on is due to the mo9on along the z-axis, which is not accounted for in our model, leading to a measured step length shorter than the actual one ( Fig. 3B, le9).
The overes)ma)on is likely due to par/cles bouncing against the boundary at z coordinates where the spherocylinder (x,y) sec&on has a smaller perimeter (Fig. 3B, right). (iii) Finally, when approxima&ng the solu%on in a bounded domain one must compute all trajectories that lead from the ini.al posi.on p 0 to the final posi-on p f a"er 0, 1, 2, . . ., N bounces. Compu-ng all these trajectories analy-cally is cumbersome, and therefore we limited our analysis to the 1-bounce case. This can be limi0ng in the case of fast diffusion in small compartments: when the square root of the mean square displacement is much higher than the size of the compartment, or when the acquisi5on 5me is very long, a par5cle could bounce against the surface mul2ple 2mes over the acquisi2on period.
Despite these caveats, we have shown and verified that the proposed ``folding approach'' is reasonable and, in fact, the shortcoming iden1fied above are merely computa1onal: (i) Complicated geometries can be approximated by the union of convex sets, similar to a triangula0on of smooth objects. For each convex set one can poten,ally tune and adapt our ``folding'' approach. (ii) In three dimensions one can extend the descrip-on for the rectangular plane to diffusion in cubes. Consequently, if the three dimensional confined geometry can be approximated by triangulated cubes, one can s6ll apply our ``folding'' argument. (iii) Analy6c computa6on of trajectories is computa6onally expensive. However, we have shown that based on a few bounces, the es&ma&on of the diffusion coefficient already improves significantly. If one pairs this with an approxima7on of complicated geometries by simpler ones, one can s"ll make sufficiently good approxima"ons.

A simulation-based solution to the limitations of confined diffusion
We then reasoned that an approach based on diffusion simula4ons could have advantages compared to mathema&cal models, as the shape of the compartment or the number of bounces against the surface would not be limi-ng for the outcome. In brief, we developed a method to recursively es1mate the diffusion coefficient with an algorithm that makes use of Smoldyn 46 and the SMdM technique 16,35 (Fig.  4A). Smoldyn allows the simula"on of the mo"on of par"cles within a compartment. The compartment can be either mathema,cally described, or an input of triangulated coordinates of any desired shape. Confinement is accounted for in the mo2on of the par2cles, which reflect off of the boundaries without losing velocity. We generated diffusion simula5ons with Smoldyn inside a spherocylinder, and as an#cipated the apparent diffusion is underes#mated, especially in regions close to the boundaries ( Fig.  2A, 2E, 3C). The extent of the underes'ma'on is propor'onal to the diffusion coefficient, and inversely propor$onal to the size of the confined space and to the acquisi$on $me. Our algorithm yields a diffusion coefficient that is homogenous throughout the whole compartment and matches the predefined value used to create the diffusion simula4ons.
The main steps for the opera/on of the algorithm are: (i) For simula/ons, an input diffusion coefficient (D input 0 ) is used to simulate the diffusion of par3cles in a spherocylinder, as previously described 16 . For in vivo datasets, diffusing par/cles are measured via stroboscopic illumina/on microscopy, as previously described. The diffusion map is experimentally obtained by SMdM 16,35 ; (ii) The diffusion map yields the total number of displacements per cell and the measured diffusion coefficient (D output 0 ) for every posi%on; (iii) The star%ng (x,y) coordinates of all the observed par3cles are used to place them inside a simulated spherocylinder, and their z coordinates are randomly assigned. This is done by taking a value from a uniform random distribu.on ranging from the lowest to the highest z value that a par+cle can have at that specific (x,y) posi'on inside the spherocylinder; (iv) The star'ng posi'ons (x,y,z) of the par$cles that belong to a specific pixel of the original diffusion map are selected; (v) Diffusion simula'ons are run, using as input diffusion coefficient the D output 0 value obtained via SMdM; (vi) The diffusion observed by SMdM is analyzed to obtain a new diffusion coefficient for the specific pixel (D output 1 ); (vii) The absolute difference (ε) between D output 0 and D output 1 is calculated; (viii) The diffusion simula'on process (and subsequent SMdM analysis) is repeated recursively, every 'me using a different The difference (ε) between D output 0 and the one obtained through The recursion is run un4l the difference between D output N and D output 0 reaches a minimum; (xi) The input value (D input N-1 ) that led to the minimal difference between D output N and D output 0 is then used to create a new map, which carries D input N-1 in the posi*on of the analyzed pixel; (xii) The process is repeated for every pixel, un9l a map is obtained in which the diffusion coefficients represents the unbiased diffusion values, that is diffusion values that are corrected for the effect of confinement and for the underes1ma1on caused by a two-dimensional representa.on of a three-dimensional mo+on. (xiii) The whole rou+ne is repeated ten +mes for each cell. The output maps obtained in each itera-on are averaged to account for the randomness introduced when assigning the z star%ng posi%on to all par%cles and for the randomness introduced by Smoldyn in simula%ng the diffusion of par-cles 46 . (Fig. 4A).

Simulation-based Reconstructed Diffusion overcomes limitations of confinement caused by cell boundaries
Simula'on-based Reconstructed Diffusion (SbRD) together with SMdM allows obtaining more accurate diffusion maps. It is possible to retrieve the actual diffusion coefficient also for the regions close to the cell boundaries (Fig. 4B). We also simulated scenarios of a cell displaying slower diffusion at one cell pole, and observed that SbRD is correctly iden6fying the region with slower diffusion and higher diffusion coefficients in the rest of the cell ( Supplementary Fig. 2).
We then benchmarked SbRD against SMdM by varying D input from 0.01 to 110 μm 2 /s and analyzing D output in the innermost 100 nm 2 square region of the simulated spherocylinder, which is the region least affected by the effect of confinement, and in a cell pole, which is the region most affected by the effect of confinement. We observe for the innermost region that D output obtained via SMdM decreases to 90% of its input value already for D input of 10 μm 2 /s, and that the devia-on increases as expected with D input (Fig. 4C, top panel). The decrease is more pronounced when the cell pole is analyzed (Fig. 4C, bo>om panel), with D output obtained via SMdM decreasing to 90% of its input value already for D input of 1 μm 2 /s. Importantly, the output of SbRD remains stable throughout the whole set of measurements, with diffusion coefficients near 100% of D input (Fig. 4C). Billiard fitting of rod-shaped bacteria The shape of E. coli cells is generally assumed to be a spherocylinder [48][49][50][51][52] . So%ware is available to determine the shape of E. coli cells from two-dimensional images 53-55 , but es(ma(on and subsequent triangulariza%on of their three-dimensional shape does not necessarily lead to the true shape of the cell. For instance, invagina/ons or protuberances observed on the xy plane might not be observed on the xz plane. By analyzing our microscopy data, we observed that the vast majority of cells had a shape that could be very well approximated by a two-dimensional projec/on of a spherocylinder. Therefore, we decided to use this shape for cell detec0on and modeling, which at the same 0me allowed for an easy implementa#on in Smoldyn. To apply SbRD to microscopy data, we developed an algorithm to automa&cally detect cells as billiards in microscopy images (Fig. 5A). Firstly, each field of view was filtered for background noise, yielding clusters of points represen6ng cells. Point clouds were rotated so that their major axis was aligned to the x-axis, and subsequently clustered using the equa4on of a billiard (Eq. 8).
We refine the cell selec,ons via Maximum Likelihood Es,ma,on using the following assump,ons (see Methods -cell clustering and detec.on): (i) fluorescent points are uniformly distributed throughout the cell; (ii) the cells, here E. coli, are modeled as spherocylinders; therefore the 2D projec7on of their fluorescence (shape of a billiard) appears more populated in the center than near the cell boundary; (iii) all the fluorescent spots observed in a cell have the same probability of being noise; and (iv) every fluorescent point is equally likely to be noise or to be a photoconverted mEos3.2. We then filtered the selected cells for a minimal length of 0.65 µm and a maximal width of 1.5 µm (see methods -cell clustering and detec.on). If the length of the final billiard describing the shape of the cell was bigger than 3 µm, the algorithm separates the cluster into two billiards, each having half the length of the original one. The refinement step was then repeated for every cluster of two-billiards. This allowed the accurate detec)on of newly divided cells (Fig. 5B), which is o<en not possible with standard clustering methods, such as Voronoi (Fig. 5C). In our previous work 16 we acquired a dataset of diffusing proteins of different complex mass (Supplementary Table 1) using SMdM. Here, each cell was clustered using Voronoi clustering. We re-analyzed the full dataset from our previous work with the new clustering method. Importantly, we observe no significant difference between the two datasets, both in the number of detected cells and in the diffusion coefficients obtained via SMdM analyses for the cell center (Fig. 5D, top) and for the cell poles (Fig. 5D, bo9om). We then used the informa"on of each cluster to recreate a spherocylinder in Smoldyn 46 having the shape of the corresponding cell, to which we applied the SbRD algorithm (see "A simula5on-based solu*on to the limita*on of confined diffusion" from point iii to point xi) to reconstruct diffusion maps, corrected for the confinement effect and bias by 2D models to describe a 3D mo/on, of the single-molecule fluorescent microscopy data. SbRD correlates the confinement-corrected diffusion coefficient with the perceived viscosity of the cytoplasm An advantage of our clustering method is the possibility to precisely iden7fy the cell poles and the cell center by using the radius of the billiard (Fig. 5B), allowing us to analyze these regions separately for every cell. We then compared the results obtained via SbRD with the results obtained via SMdM on our previous dataset 16 (Supplementary Table 1), and using the new clustering method for cell detec<on ( Table 1, Fig. 6). We observe a significant difference in the observed diffusion values for faster diffusing proteins, while for the slower diffusing par4cles the difference appears to be not significant. Notably, we observe a higher sta-s-cal significance in the difference between the diffusion coefficient observed at cell poles. These results are in line with the observa4ons made via simula4ons, which indicate that (i) given the fast acquisi/on /me used in SMdM, the effect of confinement on the measured diffusion results less pronounced for slower diffusing par2cles; and (ii) that the confinement effect is more pronounced in the cell poles than in the cell center.
We find a correla,on between the diffusion coefficient and the complex mass, while we do not observe any correla*on between the diffusion coefficient and the number of interac4ons of the analyzed proteins ( Supplementary Fig. 3). However, since SbRD yields diffusion values corrected for confinement effects, we obtain a new and improved correla4on between diffusion and complex mass, with the diffusion coefficient scaling as D = αM -0.6 . Consequently, the proposed correla4on between diffusion and perceived viscosity (η) 16 changes to η = αM complex 0.27 (Supplementary Fig. 3).

Fig. 6. Comparison of the diffusion values obtained via SMdm and via SbRD.
Comparison of the apparent diffusion coefficient obtained via SMdM and of the confinement-corrected diffusion coefficient for both the cell center (top) and the cell poles (bo1om) for the dataset of proteins tagged with mEos3.2 16 . Asterisks indicate sta*s*cal significance obtained via a Mann-Whitney U test for non-normally distributed data.

SbRD versus SMdM analysis of the diffusion coefficient at the cell poles
Diffusion measured near the cell boundary and in the cell pole regions of rod-shaped bacteria appears slower than in the cell center due to confinement effects 16,38,[56][57][58] . We recently showed 16 that the ra'o between the diffusion coefficient at the cell poles and cell center is lower for SMdM data than for simulated data, where par0cles are treated as mathema0cal points moving of random mo0on. This indicates that the slowdown observed in cells must be due to some physiological effects, such as increased crowding in the polar region, possibly due to aggrega5on of old or damaged proteins, the presence of the transla.on machinery, or dynamic structures genera!ng over-and undercrowded regions 16 . One of the key unanswered ques3ons about the diffusion measured at the cell poles is: how much of the observed slowdown is due to confinement and how much is due to physiological effects? It is not possible to decouple these effects by SMdM or other single molecule microscopy techniques. The ra#o between the diffusion coefficient at the cell poles and the cell center obtained by SbRD correlates linearly with the ra,o obtained by SMdM (Supplementary Table 1). This is not observed in simulated cells (Fig. 7A). Therefore, the slowdown observed at the poles cannot be a<ributed solely to the effect of confinement. We compared the ra0o between the diffusion at the cell poles and the diffusion at the cell center obtained via SbRD with the one obtained via SMdM for all the analyzed cells clustered as billiards.
We observe a D pole /D center ra#o of 0.74 ± 0.13 for SMdM and 0.80 ± 0.16 for SbRD (Fig. 7B). We analyzed the difference with a Mann-Whitney U rank test for non-normally distributed data and obtained a pvalue << 0.01. We therefore conclude that about 20% of the previously observed 25% slowdown at the cell poles can be a+ributed to confinement effects.

Effect of Rifampicin and Erythromycin on diffusion at the cell poles
The slowdown in diffusion at the cell poles observed via SbRD can be due to: (i) the presence of polysomes, large structures composed of several ribosomes bound to the same mRNA molecule; (ii) clusters of aggregated or damaged proteins that hinder the mobility of other molecules; (iii) or the presence of unknown substructures. To test the first hypothesis, we treated E. coli cells expressing mEos3.2 with the an0bio0cs rifampicin and erythromycin, which disrupt the polysomes via different mode of ac(ons. Rifampicin inhibits DNA-dependent RNA biosynthesis by inhibi0ng the bacterial RNApolymerase 59,60 , which leads to rapid RNA deple3on, par3cularly of mRNA 61,62 , while erythromycin inhibits the assembly of the large ribosomal subunits 50S 63 . Hence, the use of these an/bio/cs should make the diffusion coefficients at the poles and middle of the cell similar if the polysomes would form a major hindrance for the mobility of mEos3.2. However, we do not observe a significant change in ra=o of the diffusion coefficients (Fig. 7C), with values of 0.84 ± 0.12, 0.86 ± 0.12 and 0.82 ± 0.12 for untreated cells, cells treated with erythromycin and cells treated with rifampicin, respec4vely. This is confirmed by the Mann-Whitney U rank test for non-normally distributed data. To further inves5gate the effect of the an#bio#cs, we analyzed the absolute values of diffusion. For cells treated with erythromycin, the diffusion at the cell center and at the cell poles show a moderate increase compared to the control sample, while for cells treated with rifampicin the increase is much larger (Fig. 7D)

Analysis of diffusion at the cell poles indicates asymmetry that correlates with aging
We then analyzed the microscopy data for differences between the cell poles ( Supplementary Fig. 4). We reasoned that differences in aging of the two poles could lead to differences in diffusion, especially because misfolded and aggregated proteins tend to accumulate at the old pole [19][20][21] . We acquired SMdM data of newly divided cells by selec1ng fields of view with cells that had just completed the division process (Fig. 5B, 5C), and we subsequently applied SbRD to the datasets. We determined the diffusion coefficient corrected for confinement effects of mEos3.2 in each cell individually and find that the diffusion at the old pole is significantly slower than at the new pole. The ra7o between the diffusion at the new cell pole and the cell center is 0.86 ± 0.15, while D old-pole /D center is 0.80 ± 0.13 ( Supplementary Fig.  5). With the assump/on that the new cell pole enables on average faster diffusion than the old cell pole, we subtracted D old-pole from D new-pole for each cell. We obtained a distribu3on of residual diffusion coefficients with a mean higher than zero (Fig. 7E). We performed a one-sided Wilcoxon signed-rank test to confirm whether the observed difference was significant, and obtained a p-value < 0.05. These findings confirm the hypothesis that aging influences the structure of the cytoplasm at the poles of E. coli, causing macromolecules at the older cell pole to diffuse slower than at the new one.

Discussion
We developed a new method to obtain diffusion coefficients that are not affected by confinement effects and bias by 2D modeling of a 3D mo2on. This method is key not only for measuring lateral diffusion in small compartments, but also for diffusion in proximity of boundaries, such as the plasma or organellar membranes. Despite the enormous advancements offered by SMdM in probing molecule mo:on compared to other methods, the values obtained near boundaries are affected by the effect of confinement. Our newly developed method SbRD allows examining more precisely regions of the cell that are small and or geometrically more complex.
A method to reconstruct confinement-corrected diffusion coefficients in small cells The main advantage of the SbRD method is that the shape of the analyzed compartment does not limit the analysis. In fact, compartments of any shape, that can be visualized and reconstructed via triangulariza*on, can be used to recreate an iden*cal virtual compartment. In this way, our method allows reconstruc-ng unbiased diffusion in heterogeneous compartments such as those in eukaryo-c cells.
We show by simula/ons that lateral diffusion measurements performed in small compartments, such as the prokaryo+c cell, are bound to underes+mate the diffusion coefficient, not only in the regions near the boundary but also in the cell center. The confinement effect in the cell center is mostly due to 2D observa(on of a mo(on in 3D. Using the here developed tool for cell clustering, we precisely detect spa$al informa$on such as radius, length, center and orienta$on angle of a cell. We use this informa$on to reconstruct cells with an assumed spherocylindrical shape in Smoldyn 46 . Recursive diffusion simula!ons then yield the diffusion coefficient corrected for the confinement effect and the spa!al component of mo*on. This approach led us to (more precisely) es*mate the dependence of the diffusion coefficient on the complex mass of the diffusing species and infer from the data the viscosity perceived by a molecule of given molecular mass.
We observed a linear dependence of the ra3o between the diffusion at the cell poles and the cell center for data acquired via SbRD and via SMdM, as shown in Fig. 7A, indica%ng that the observed slowdown in diffusion at the cell poles can be a1ributed to physiological effects, and not solely to confinement.

Asymmetric diffusion at the cell poles correlates with aging
We previously obtained indica3ons that the diffusion at the cell poles is slower than in the center of E. coli cells 16 . We now determine precisely how much slower the diffusion is at the cell poles, and we show that disassembly of polysomes and deple1on of mRNA by an1bio1c treatment do not affect the differences in diffusion between the poles and the center of cell. In fact, the diffusion coefficient increases by a similar percentage in each region of the cell, sugges5ng that the an5bio5c treatment has decreased the overall viscosity of the cell.
Moreover, we were able to precisely detect dividing cells and the new cell poles, and we show that the ones at the division site exhibit a faster diffusion compared to the old cell poles. The diffusion coefficients at the new and old pole are 86% and 80%, respec-vely, of the value measured at the cell center. In eukaryo.c cells aging is accompanied by an increased cytosolic crowding 64 . Here, we hypothesize that rela.ve slowdown of diffusion at the old pole is consistent with an increase in crowding and an indica'on of aging in E. coli.
Previous studies suggested the possibility of accumula6on of aggregated proteins at the cell poles of E. coli as a possible cause for the observed aging effects 19 . In the recent study from Łapińska et al. 45 , however, the authors do not observe protein aggrega3on in the form of inclusion bodies and they argue that proper techniques for the inves,ga,on of the structure of cell poles are lacking. Here, we developed a method to more precisely assess the structure of the cytoplasm in the cell pole regions. We also do not see indica)ons of protein aggrega)on such as inclusion bodies in our dataset. We conclude that the slower diffusion observed at the old cell pole is an indica4on of the presence of par4ally aggregated misfolded macromolecules.
The E. coli cell cycle is divided in three different phases: the B period, which goes from the cell birth to the beginning of DNA replica3on; the C period, which goes from the beginning to the termina3on of DNA replica,on; and the D period, which represents the ,me between the termina,on of DNA replica(on and cell division 65,66 . During the D period, the two chromosomes segregate to two different parts of the cell, and a period of protein synthesis is necessary for comple3ng the cell division process. If accumula&on of sta&c or semi-sta$c structure is a sole characteris$c of the cell poles, then the newly formed cell pole should exhibit a diffusion value close to the one measured at the cell center. Since we observe a slower diffusion at the new cell pole (albeit less slow than at the old pole) compared to the one measured at the cell center, we conclude that the accumula1on of damaged or aggregated protein is not a specific characteris.c of the cell pole, rather it is a consequence of steric hindrance exerted by the nucleoid, which causes accumula.on of structures in nucleoid-free regions of the cell. Based on the currently available literature 19,44,45,67 , we hypothesize that the observed slowdown is maintained throughout the cell cycle, and that it is possibly passed down from mother cell to daughter cells, where one of the daughter cells will inherit the "slow pole" from the mother, while the other will inherit the "fast pole" (Fig. 8). As the cell cycle progresses this difference is maintained. When the cell divides and the septa5on ring forms (red), the two daughter cells will inherit the two cell poles. One of the cells will have the oldest pole as its old pole, while the other will have the new pole of the mother as its old pole.

Concluding remarks
We developed Simula.on-based Reconstructed Diffusion (SbRD) to determine diffusion coefficients in any compartment, corrected for confinement effects and by the mo2on of par2cles along the z-axis. We applied the technique to a previously recorded dataset. Using SbRD, we obtain a more precise correla'on between the diffusion coefficients and the complex mass of the diffusing species, from which we infer the perceived viscosity for proteins diffusing in the cytoplasm. We recorded new singlemolecule displacement datasets to characterize the slower diffusion at the cell poles, and to determine differences in confinement-corrected diffusion at the old and new pole. We correlate slower diffusion at the old poles with aging of the cells. We argue that our method and analyses provide new possibili9es for inves*ga*ng the mechanism of aging of bacteria and other types of cells.

Live-cell single-molecule microscopy
Media prepara)on, cell culturing, measurements setup and live-cell imaging was performed as described 16 . Briefly, for each experiment we started a pre-culture of E. coli, bearing a pBAD plasmid for the expression of mEos3.2, by scratching a glycerol stock with a sterile inocula<on loop and dipping it in a 14-mL plas(c culturing tube containing 3 mL of LB medium, prepared following the formula of 10/10/5% (w/v) in MilliQ of NaCl, tryptone (Formedium), and peptone (Formedium), respec3vely, and supplemented with ampicillin (100 µg/mL). We incubated the pre-culture overnight at 30°C, with shaking at 200 rpm. On the following day we transferred 30 µL of the LB pre-culture into 3 mL of MOPS-buffered minimal medium (MBM), prepared following the formula in 68 , supplemented with 0.1% (v/v) glycerol and ampicillin (100 µg/mL). Cultures were incubated overnight at 30°C, with shaking at 200 rpm. The next day cells were diluted to a final OD 600 of 0.05 to 0.08 into prewarmed MBM containing 0.1% (v/v) glycerol, ampicillin (100 µg/mL) and 0.1% (w/v) L-arabinose, and incubated the at 30°C, with shaking at 200 rpm for 4 to 6 hours before microscopy experiments. Right before the measurements, the cultures were spun down in a tabletop centrifuge and concentrated three 4mes in the growth medium.
To ensure a constant temperature of the microscope during the imaging process, the instrumenta3on was turned on 4 to 5 hours before the measurement, to minimize the xy dri$ of the samples. Cells were imaged on a clean, non-func%onalized high-precision glass slide {specs, manufacturer}, previously sonicated in 5M KOH for 45 minutes and then rinsed 10 8mes with MilliQ, followed by a drying process via pressurized air. Immobiliza2on of the samples was achieved by deposi2ng 5 µL of concentrated cell suspension on the glass slied and then pressing the cells against the glass surface with solidified agarose pads having the same composi0on of the MBM medium with a final concentra0on of agarose of 0.75% (w/v), formed inside a polydimethylsiloxane (PDMS) chamber.
Once the cells se*led, we selected an area of our field of view to perform the measurements. We adjusted the focus and the laser beam angle to obtain the highest number of foci, which resulted in the beam angle slightly below that of the cri3cal angle for total internal reflec-on [highly inclined and laminated op,cal sheet microscopy 69 ]. The camera and the laser were then synchronized in the stroboscopy mode, with illumina5on pulses necessary to first photoconvert and then detect mEos3.2 every 1.5 ms 16,35 . For a detailed overview of the scripts used for managing the microscope, we refer to our code 70 .
Erythromycin treatment was performed for 1 hour a5er cells reached an OD 600 of 0.12-0.15 on the day of measurement, by adding erythromycin to the cell culture to a final concentra5on of 250 ng/µL. Agarose pads were supplemented with the same erythromycin concentra3on.
Rifampicin treatment was performed as described for erythromycin, using a final concentra'on of rifampicin of 500 ng/µL. Agarose pads were not supplemented with rifampicin, as this influenced the photoconversion of mEos3.2.
To analyze dividing cells, we visually inspected different fields of view in every sample, and selected areas in which a pair of obviously dividing cells were observed. Dividing cells were then analyzed separately during data analysis.

Cell clustering and detection
The first step in our data analysis pipeline is represented by the detec4on and clustering of the cells in the imaged field of view. First, the whole analyzed field of view was converted into a 2D histogram in which every bin represents a certain number of fluorescent points. The background value was calculated by taking the median value of the bins. Only the fluorescent points belonging to bins having a value higher than the background value plus one standard devia4on were kept and used to create point clouds. For each point cloud the eigenvectors were obtained through the calcula6on of the covariance matrix, which allowed calcula0ng the angle between the first eigenvector and the x axis. From this, an appropriate rota)onal matrix was applied to the xy coordinates of the point cloud, to align its major axis parallel to the x-axis. This allowed obtaining a first set of features of the point cloud, namely the length, the width, the rota,on angle and the center, which we used to describe a billiard encompassing the point cloud. Shape refinement was obtained by fi7ng an improved billiard around the point cloud via maximum likelihood es.ma.on. The fi3ng was performed by applying the following assump.ons: (i) the fluorescent molecules are uniformly distributed throughout the cells. Since the cells are spherocylinders imaged in two dimensions, the number of molecules observed is directly propor2onal to the thickness of the cells; therefore the boundary areas are less populated than the center; (ii) every observed point is equally likely as any other to be due to random noise; (iii) the probability of a point being random noise is equal to the probability of a point being a fluorescent molecule. From these three assump8ons we obtain the following probability mass func5on for each par5cle (Eq. 9): Where N is the number of par0cles (X,Y) inside the spherocylinder, V is the volume of the spherocylinder, which is modeled based on all the points (X,Y) within the spherocylinder, and h x,y is the thickness of the spherocylinder at the xy coordinate of the detected point. All these parameters depend on the size of the spherocylinder, which is described by its length, its radius, its center and its rota5on angle. Therefore, these are used as fi/ng parameters to iden5fy the best spherocylinder describing the detected point cloud.
The iden(fied cells are then filtered based on their shape, discarding cells that are shorter than 0.65 µm (possibly cells that are par0ally out of the field of view), or wider than 1.5 µm (possibly noise or dri#). Cells having a length bigger than 3 µm were automa5cally reanalyzed as dividing cells. In case of overlapping billiards for dividing cells, the intersec4on points were iden4fied by calcula4ng the intersec(on of the two semicircles describing the two adjacent cell poles. From this, it was possible to calculate the volume of the two spherical sec2ons of the neighboring cell poles. Since any of the observed points could belong to either cell, the final volume used to calculate the probability density had to be adjusted by adding the intersec2on volume (Eq. 10): Where N is the number of par0cles (X,Y) inside the two spherocylinders, V(X 1 ,Y 1 ) is the volume of the first spherocylinder, V(X 2 ,Y 2 ) is the volume of the second spherocylinder, V(X int ,Y int ) is the volume of the intersec(on between the two spherocylinders, h 1 is the thickness of the first spherocylinder at the xy coordinate of the detected point, and h 2 is the thickness of the second spherocylinder at the xy coordinate of the detected point.
Finally, the iden-fied cells are visually inspected and discarded if not suitable for analysis (e.g. dividing cells not correctly iden-fied, for which the total length is shorter than 3 µm). The fi6ed spherocylinders are then used to create clusters from the point clouds. The data analysis is then performed on each cluster separately, ignoring the points that are not included in the cluster. More informa(on on the cell detec%on and clustering can be found on our code 70 .

SMdM analysis
SMdM analysis was performed as described previously 16 , with the excep+on that cell clustering was performed prior to peak pairing. Briefly, we recorded several consecu:ve movies for each field of view and paired the observed localiza2ons from the two consecu2ve frames of the stroboscopic illumina2on pa#ern. For single-molecule analysis, we used the STORM-analysis package developed by the Zhuang laboratory, which is included in the 3D-DAOSTORM program for peak detec0on 71 . A$er a full movie was analyzed, the localiza/ons were corrected for xy dri$.
All the detected peaks in each field of view were used for clustering and for finding the shape of the spherocylinder that best describes the shape of the cell, as reported in the sec3on above.
Displacements were obtained from all the peaks belonging to a single cell, by pairing localiza9ons from the two consecu+ve frames of the stroboscopic illumina'on pa+ern. We set a maximum distance of 600 nm between any two peaks to be paired: the distance is then used to find all possible peak pairs for each couple of frames. To obtain a displacement, we match each peak in the first frame of the couple with all the peaks falling within a radius of 600 nm in the second frame. This procedure is repeated for all frame couples of each field of view. A hard filter based on the number of detected displacements was then applied, discarding cells with less than 2000 or more than 20000 displacements, as described 16 .
A pixel map with pixel size of 100 nm 2 was obtained for each cell, with every pixel containing the informa(on of all the peak pairs for which the star(ng posi(on is located inside the pixel itself. Each pixel of the map containing a minimum number of displacements (set to 10 in our study) was then fi;ed using a modified two-dimensional probability density func3on (PDF), which accounts for a linear background effect k, which can be caused by an ambiguity in the assignment of peak pairs 16 (Eq. 11): Since t is known, as it represents the /me between two stroboscopic laser pulses, the PDF was fi!ed on the detected displacements r, using the diffusion coefficient D and the background value k as fi%ng parameters. Displacements were detected using the MLE clustering method (see sec9on above) for SbRD. Displacements were detected both using the MLE clustering method and Voronoi clustering when comparing the two clustering methods using SMdM.
Given our advanced cell detec/on method, we could perform an accurate iden/fica/on of the different cell regions, namely the cell poles and the cell center. Fi&ng a billiard to detect the shape of the cells allows obtaining precise informa2on about their length and their radius, which in turn allow iden2fying the cell pole regions and the central region of the cell. All the displacements belonging to the same region were then used to perform the fi2ng using equa5on 10, yielding informa5on about the diffusion coefficient in the different regions of the cell 16 .
The dependence of the diffusion coefficient on the complex mass was fi7ed using a power law rela%onship D = αM complex β , where M complex is the complex mass and α and β are fi&ng parameters 16 .
Fi#ng was performed using the func5on curve_fit included in the SciPy library 72 .

Smoldyn simulations
Simula'ons were performed using the so5ware Smoldyn 46 , as described 16 . A diffusion coefficient and a !me-step length are used as input for the simula2ons, together with the total simula2on 2me. At every !me step, Smoldyn randomly selects a step length from a normal distribu!on having as mean the squared mean squared displacement calculated from the input diffusion coefficient, as well as a random direc&on in the xyz space for each par+cle. These values are used to simulate the mo1on of every par$cle in the system at every $me step, un$l the total simula$on $me is reached. In our simula$ons we used a 'me step of 0.1 ms and a total simula'on 'me of 2 seconds. Par'cles every 15 steps (1.5 ms) were then paired together in displacements, the results were benchmarked against the microscopy data.
We used Smoldyn to generate two separate datasets. First, we simulated the mo8on of par8cles using input diffusion coefficients ranging from 0.01 to 110 µm 2 /s in a spherocylinder having length and width of 2.25 and 0.9 µm, respec4vely, as these reflect the average cell size observed in our previous work 16 .
We then generated a second dataset using input diffusion coefficients ranging from 1 to 20 µm 2 /s in a spherocylinder with a length ranging from 1.4 to 2.9 µm and width ranging from 0.6 to 1.5 µm, always keeping the ra,o between length and width higher than 2 and lower than 4 to reproduce actual dimensions of E. coli.
We analyzed the output of our simula/on using the same approach adopted for SMdM measurements and compared the results with those obtained by SbRD (see sec8on SbRD analysis). The equa*on used for simulated data does not account for linear background correc3on (Eq. 12): . 12 More informa+on about the simula+ons can be found on our code 70 .

SbRD analysis
Simula'on-based Reconstructed Diffusion (SbRD) was applied by using Smoldyn simula:ons on SMdM analyses. A collec,on of points represen,ng a cell, either coming from a Smoldyn simula0on or from a microscopy experiment, is analyzed using SMdM (see sec8on SMdM analysis). This analysis provides a pixel map of the cell, in which for each pixel the diffusion coefficient, as well as the xy start and end posi%ons of each displacement, are known. For microscopy measurements, the new cell clustering method allowed us to determine the length and radius of the spherocylinder in which the diffusing par$cles are confined, which could be modeled in Smoldyn. In the case of dividing cells, the septum is modeled as a reflec,ng wall passing through the intersec,on points of the two billiards, encompassing two separate cells, and parallel to the z axis. In the case of data from Smoldyn simula6ons, the length and radius of the spherocylinder are known. These informa5on were then used to start a recursive simula'on in Smoldyn by placing a number of par'cles equal to the number of displacements in their respec&ve xy star%ng posi%on inside the pixel, with the z posi%on randomly assigned to each par2cle inside the spherocylinder, which in the case of microscopy data was modeled as described in Cell clustering and detec.on. A simula*on las*ng for 1.5 ms, with simula*on steps of 0.1 ms is started, using as input diffusion coefficient the value of the pixel obtained via SMdM (see sec7on Smoldyn simula-ons).
The output of this simula/on is then used to perform a fi4ng using equa/on 11, with D as fi%ng parameter. The squared difference between the output diffusion obtained via simula,ons and the diffusion obtained via SMdM is then calculated. The program then recursively iterates the simula:on process un*l such squared difference reaches a minimum. The input diffusion coefficient used to obtain the output diffusion coefficient that minimizes the squared difference is then regarded as the real diffusion coefficient of the pixel. This process is then repeated 10 8mes for each pixel to account for the randomness introduced by Smoldyn 46 in the choice of the step length and the direc/on of mo/on, as well as for the randomness introduced in the placing the par4cle along the z-axis. The process is then repeated for every pixel of the original SMdM map, from which an SbRD map is obtained. The pixel-bypixel differences between the SbRD map and the SMdM map are used to construct a difference map. More informa+on about the SbRD analysis can be found on our code 70 .

Statistical analysis
All sta's'cal analyses were performed using the Python package stats from the SciPy library 72 . Shapiro-Wilk test for normality 73 was used to check whether the data are normally distributed, using a level of confidence of 5%. The test assumes the null hypothesis for data that are normally distributed. Therefore, if the obtained p-value is lower than 0.05 the null hypothesis is rejected and the data are assumed to be non-normally distributed. Non-normally distributed datasets are visualized via kernel density es3ma3on.
The Mann-Whitney U rank test 74 was used to test whether the means of two non-normally distributed datasets are equal. In the case of comparing means of datasets, i.e. when no prior assump.ons are made and no precise outcome is expected, such as in the case of comparing the diffusion coefficients of cells treated with an.bio.cs, a two sided test was performed. In the case of comparing means of datasets in which a specific outcome was expected, such as in the case of comparing the diffusion coefficient of the two different cell poles, a one sided test was performed.
The Wilcoxon signed-rank test 75 was used to test whether the median of a dataset coming from paired measurements is significantly different from zero. This test was used to assess whether the difference in diffusion between the new cell pole and the old cell pole was significantly higher than zero, therefore it was conducted as a one sided test. Sta/s/cal significance in pictures are indicated with 1 asterisk (*) for p-value < 0.05, 2 asterisks (**) for p-value < 0.01 and 3 asterisks (***) for p-value < 0.001.