extfieldsofstudy:
Chemistry
Computer Science
Physics

text:
Dynamic energy budget approach to evaluate antibiotic effects on biofilms

Quantifying the action of antibiotics on biofilms is essential to devise therapies against chronic infections. Biofilms are bacterial communities attached to moist surfaces, sheltered from external aggressions by a polymeric matrix. Coupling a dynamic energy budget based description of cell metabolism to surrounding concentration fields, we are able to approximate survival curves measured for different antibiotics. We reproduce numerically stratified distributions of cell types within the biofilm and introduce ways to incorporate different resistance mechanisms. Qualitative predictions follow that are in agreement with experimental observations, such as higher survival rates of cells close to the substratum when employing antibiotics targeting active cells or enhanced polymer production when antibiotics are administered. The current computational model enables validation and hypothesis testing when developing therapies.

Introduction
Biofilms are bacterial aggregates that grow on moist surfaces, encased in a self-produced polymeric matrix, see Figure 1. The matrix creates a favorable environment for their development, facilitating nutrient, oxygen and waste transport [30]. It also acts as a shield against external aggressions by flows, disinfectants and antibiotics. The minimal bactericidal concentration (MBC) and minimal inhibitory concentration (MIC) of antibiotics to bacteria in their biofilm habitat may be up to 100-1000 fold higher compared with planktonic bacteria [1,16].
Implant associated infections typically involve biofilm growth on the surface of the implant [40]. They form on medical equipment and prostheses, such as pacemakers and endotracheal tubes, central lines, intravenous catheters, stents and artificial joints. Bloodstream infections, and many other hospital-acquired infections, may be caused by them. Biofilms may also spread on body surfaces such as heart valves (endocarditis), teeth, the lungs of cystic fibrosis patients (pneumonia), the middle ear and nose (otitis, rhinosinusitis), bones (osteomyelitis) or in chronic wounds [16]. The biofilm matrix hinders phagocytosis and other actions of the inmune system. Bacteria surviving standard antimicrobial therapies are able to reproduce, originating chronic infections [38]. To tackle this problem, we must be able to understand how antibiotics act and how resistance to antibiotics develops.
Antibiotics affect cells in diverse ways [19]. β-lactams (penicillins, cephalosporins, carbapenems) and glycopeptides (vancomycin) inhibit cell-wall synthesis. Aminoglycosides (streptomycin, gentamycin, tobramycin) inhibit protein synthesis. Quinolones (ciprofloxacin, ofloxacin) inhibit DNA replication. Tetracyclines inhibit translation. Polymyxins such as colistin disrupt charge distributions in the outer cell membrane. Eventually, the damage caused to the cell produces its death. To exert their antibacterial action, antibiotics undergo a certain process. They need to penetrate the cells, remain stable and accumulate to reach inhibitory concentrations. Sometimes they have to take an active form. After detecting their target, they interact with it to exert their action. Interferences in either of these processes may result in cell resistance to the antibiotic. Resistance typically proceeds through efflux systems (the antibiotic is pumped out of the cell), chemical alterations of the antibiotic (cellular enzymes degrade it), mutations in antibiotic target molecules, and non-heritable resistance caused by environmental conditions [9,16,19,36,37]. The main resistance mechanisms for different types of antibiotics are summarized in [19].
The biofilm environment enhances bacterial resistance in a number of ways. Biofilm development is influenced by quorum sensing [8]. Through quorum sensing mechanisms, bacteria sense when a critical number of them are present in the environment. They respond by activating genes that produce exopolysaccharides [7]. The polysaccharide matrix surrounding the bacterial community delays diffusion of antibiotics inside the biofilm. Nevertheless, direct measurements suggest that some antibiotics equilibrate within the biofilm [6] after a waiting time. Pseudomonas aeruginosa tends to be the main source of gram-negative infections in intensive care units in developed countries [32]. P. Aeruginosa and other bacteria express β-lactamase, an enzyme that attacks β-lactams. An enzyme breaking the antibiotic at a rate at which it crosses the cell membrane combined with delayed diffusion might explain resistance to penicillins, but not to other antibiotics [6].
As mentioned above, as we penetrate from the outer biofilm surface towards the interface with the substratum, gradients of oxygen and nutrients develop. Oxygen depletes [5,42]. These gradients result in increased doubling times for cell division and reduced bacterial metabolic activity. The intensity of metabolic processes is stratified: high activity in the outer layers and slow growth or no growth in the inner core. These dormant cells are partially responsible for tolerance to antibiotics. Popular monotherapies with β-lactams are only active against dividing cells [3], forcing combinations with antibiotics that are active against nondividing cells, such as colistin [16]. Oxygen limitation and metabolic rates are also important factors enhancing the tolerance of biofilms to ciprofloxacin and aminoglycosides [41].
In the biofilm, bacteria are exposed to oxidative stress, that causes hipermutability. Enhanced production of reactive oxygen species (ROS), either released in response to the infection or produced by the alterations in the DNA repair system of the bacteria, leads to an environment with low oxygen tension filled with oxygen radicals [16,29]. Augmented β-lactamase synthesis, overexpression of efflux-pumps and increased EPS production follow [4,16,26]. Studies with toxicants have also shown increased EPS production and ability to adapt to the toxicant, repairing damage to the cell [13]. Quorum sensing inhibitors [15], efflux inhibitors [27,28], antioxidants reducing the oxidative stress and mutations [29], together with enzymes able to dissolve the biofilm matrix [2], may provide strategies to overcome resistance mechanisms.
Mathematical modeling can assist in the design of therapies and the interpretation of experimental data. Early models were able to reproduce elementary qualitative behavior. Reference [35] uses coupled reaction-diffusion equations for the concentration of oxygen, antibiotics and the volume fractions of live and dead cells to predict survival profiles inside thick biofilms due to slow growth. Reference [33] predicts that the biofilm matrix can not prevent diffusion of β-lactam antibiotics into the bacteria provided the amount of chromosomal β-lactamase is low. The diversity of the mechanisms involved in biofilms resistance to antibiotics suggests the opportunity of adapting dynamic energy budget (DEB) frameworks to describe the effect of antibiotics on them. DEB models have already been exploited to describe the effects of toxicant exposure on populations of floating bacteria [21]. Like antibiotics, toxicants interfere with the metabolism of cells and increase the energy required for cell maintenance. The cell requires additional energy to expel the toxicant and repair damage caused by toxicant activity (DNA, RNA, protein repair).
Dynamic energy budget models relate biomolecular processes to individual physiology and population dynamics [23]. They have been successfully used to study scaling behaviors of all sorts of living beings, from plants and animals to cells [14]. Basic DEB models describe acquisition of biomass and energy, as well as energy allocation for cell maintenance, growth and division. A standard description of aging and death processes is also available [23]. Surplus reactive oxygen species (ROS) causes irreparable damage to the cell. Damage components become 'damage inducing compounds', resulting in cell malfunctioning and further damage, raising the death probability. Death and damage are represented through 'hazard' and 'aging acceleration' variables. Antibiotic impact, like toxicant consequences, must be described modifying standard fluxes and rates. To represent the effect of toxicants on floating bacteria, Reference [21] included acclimation to the toxicant effect, environmental degradation due to cell products, toxicant induced ROS production, and their influence on assimilation, growth, aging and harzard rates. The only variable states characterizing the environment are substrate and toxicant concentrations.
In this paper, we analyze antibiotic effects on Pseudomonas aureginosa biofilms by means of an adapted DEB model. Unlike typical DEB models for cells [21], we introduce here spatial variations in the concentration equations to study the local influence of the biofilm shape and composition. Spatial gradients in the concentrations of oxygen and antibiotics are relevant to the bacterial survival profile in these communities. We couple numerically the metabolism of individual cells to the diffusion equations describing concentrations and some additional field. Our model accounts for EPS production, identified as a cause of antibiotic resistance. The antibiotic acts on growth, maintenance, aging and hazard rate in a similar way to toxicants. EPS production increases acclimation, while reducing growth, damage and hazard rates, see Figure 2.
We have fitted the model for a quantitative agreement with the death rates reported for P. aureginosa under different antibiotics [17] and explored the spatial variations in the dead/alive cell distribution. Our simulations reproduce expected qualitative trends. Augmenting the dose, the number of dead cells increases. Increasing the EPS extent, the number of dead cells decreases. Antibiotic presence enhances EPS production. When we apply antibiotics targeting active cells a necrotic region progressing from the outer biofilm layers to the inner biofilm core appears. Instead, antibiotics targeting slowly growing cells would destroy first the inner dormant core.
The paper is organized as follows. Section 2 describes our spatially varying DEB model. Section 3 presents the computational setting. Parameters are fitted to measured death rates in Section 4. Section 5 discusses numerical simulations illustrating the above mentioned qualitative behaviors. Finally, Section 6 states our conclusions.

Basic DEB model for a biofilm including EPS production and antibiotics
In this section we propose the basic equations for a dynamic energy budget theory (DEB) of a biofilm including polymer (EPS) production, coupled to diffusion of oxygen and antibiotic concentrations. We also introduce equations describing the effect of the antibiotic on the bacteria. The framework is similar to that developed in [21] for the study of the effect of toxicants on a homogeneous population of floating bacteria. Differently, in a biofilm we must incorporate EPS production, as well as diffusion processes in a number of spatially varying magnitudes, taking into account the specific action mechanisms of antibiotics. Also, the population of bacteria differentiates in several types. Only some of them become EPS producers. The remaining bacteria grow and divide normally, unless resources are so scarce to trigger deactivation or damage is large enough to kill them.

The basic energies
DEB is essentially a scaling theory of different types of energies characterizing the biofilm. The basic differential equations describe the fluxes of these energies. We apply a variant of the κ rule [23], see Figure 2. A κ fraction of the energy is used for bacterial growth and division. Notice that bacteria do not mature and age, just grow and divide. A 1 − κ fraction of the energy is used for EPS production.
Each alive bacterium in the biofilm evolves according to the following system of equations, setting κ = 1 for normal cells and 0 < κ < 1 for EPS producers.
1. Scaled energy density e(t): where C o is the oxygen limiting concentration, K o the limiting concentration half saturation value, f is the scaled functional response, ν is the energy conductance and ν takes into account toxic effects on conductance through ν A , defined in equation (7).
where r is the bacterial (cell biomass) production rate, m κ is the maintenance rate and g κ is the investment ratio. The remaining parameters and magnitudes are linked to environmental toxicity, in our case antibiotic concentration, see below. h is the hazard rate, a is the acclimation energy density and a M is the target acclimation energy. The symbol + stands for positive part, which becomes zero for negative expressions. Bacteria are rod-shaped, and grow in length. If we fix an average radius, equation (2) provides an equation for the time evolution of the cell length l.

Dimensionless volume of EPS
This equation defines a rate of EPS production r e . In absence of acclimation and hazard effects (a = a M , h = 0), we set the reference rate r e = kr + k , that can be fitted to experiments, as we explain in Section 5.
The dimensionless magnitudes e, v and v e represent the fundamental energies characterizing the cell (bacterium) and its EPS production. Their fluxes are written here for each individual cell, e being the total energy available for both growth and reproduction and for EPS production. v keeps track of the volume of the cell as it grows, this is the energy available for reproduction by division. v e keeps track of the energy available to form the matrix that constitutes the biofilm and protects the bacterium. They may be related to dimensional volumes V and V e multiplying by a characteristic cell volume V m . The dynamics of these variables is coupled to the evolution of other magnitudes, governed by equations we introduce next. Table 1 collects all the variables appearing in the model. Table 2 gathers the parameters, specifying units and selected values. The discussion in Section 4 suggests that a characteristic time scale for equations (1)- (3) is set by the bacterial doubling time ∼ 1 µmax or the conductance ν, in a time scale of hours. Diffusion processes for oxygen and antibiotic concentrations, however, evolve in a time scale of seconds, as concluded from the equations presented next and the parameter values in Table 2.

Diffusion of oxygen and antibiotics
Let C o (x, t) denote the concentration of oxygen and C a (x, t) the concentration of antibiotic. These quantities are described by two diffusion equations. The antibiotic accumulated inside the cell is described using the antibiotic cellular density [C IN ](t), as discussed next.
1. Oxygen concentration diffusion inside the biofilm: where µ max is the maximum specific growth rate. k is the growth associated polymer formation rate coefficient and k is the non-growth associated polymer formation rate coefficient.   2. Antibiotic concentration diffusion inside the biofilm: where R e and R are evaluated averaging r and r e in the control volumen V T . d a is the antibiotic diffusion coefficient. Equation (5) is supplemented with boundary condition C a = C a,out at the interface with the antibiotic providing fluid and the no-flux condition ∂Ca ∂n = 0 at the biofilm/substratum interface. 3. Antibiotic cellular density: where k I A is the antibiotic influx coefficient, and k O A is the antibiotic efflux coefficient. C a should be the antibiotic concentration outside the cell, that we take equal to C a given by (5) in the control volume that contains the cell. We may set [C IN (0)] = 0 initially.
In view of the parameter values in Table 2, the chemical concentrations C o and C a will evolve much faster than the other magnitudes related to cell behavior, that is, in a time scale of seconds, not hours.
We have to complete the above equations by describing how the antibiotic and oxygen concentration act on the cell and influence its energy fluxes. We mostly follow [22,21] considering the antibiotic to be a toxicant. Besides, we introduce several new DEB variables and make small modifications to account for spatial variations.

The effect on cells of oxygen and antibiotics
The previous equations take into account the decline in the organism's capacity to acquire and use energy due to respiration and to the presence of antibiotics through a diminished conductance, reduced growth due to acclimation to the antibiotic and the hazard rate, see where K V is the noncompetitive inhibition coefficient and γ ε is the environmental degradation effect coefficient. 2. Environmental degradation ε(x, t): ν ε is the environmental degradation coefficient and ν m is the maintenance respiratory coefficient. X is the cellular structure concentration computed in the control volume V T containing N cells. R is similarly computed averaging r over this control volume. The diffusive term accounts for the fact that ε feels the spatial variations. We impose no-flux boundary conditions. 3. Acclimation energy density a(t): where a M is the target acclimation energy and + stands for the positive part. 4. Hazard rate h(t): 5. Aging acceleration q(t): where h a is the Weibull aging acceleration and s q is a multiplicative stress coefficient. 6. Aging in acclimation due to dissolved antibiotic and EPS: where k I qA is the dissolved antibiotic toxicity.

Cell division, differentiation and death
Bacteria of the Pseudomonas Aeruginosa genus are rod-like. Their volume is defined by V = πρ 2 L, ρ being the radius and L the length. Increase in volume is equivalent to increase in length, setting ρ equal to their average radius. Reference [34] reports cell lengths L between 1.67 µm and 2.75 µm for P. Aeruginosa. Diameters 2ρ vary between 0.5-0.73 µm without clear relation to the length and with a preference for 2ρ = 0.6 µm. Notice that the equation for length is formally the same as the equation for volume (2), up to a factor. We may use V m = πρ 2 L m , L m being the maximum length, to nondimensionalize.
In dimensionless variables, cell division happens provided a threshold maintenance energy e d is surpassed when a threshold length l d is reached, see Figure 4. Then the cell divides and creates daughter cells.
Nondividing cells behave in different ways according to their lengths, energies, aging and hazard rates, see Figure 4: • For intermediate energies e eps1 ≤ e ≤ e eps2 the cells produce EPS according to equation (3), provided they are not newly born. A way to keep track of their age is the aging rate q, governed by (11). It is set initially equal to zero for newborn cells. We require q ≥ q eps for the cell to have the ability to become a EPS producer.
• For low energies e ≤ e eps1 , the cell just maintains itself.
• The hazard rate, given by (10), determines when cells die. Their survival probability p(t) is governed by the equation Now with the energies e, v, v e , oxygen and antibiotic concentrations C o ,C a , [C IN ] and the variables , a, h, q, which measure how such concentrations influence the energies, we have a complete set of variables and of equations. These additional rules define the behavior of each cell.

Computational framework
For computational purposes, the biofilm is identified with a slab containing cells. We divide the slab in a grid of cubic control volumes, that we use to discretize the concentration fields evaluating their spatial changes. A control volume V T may contain several cells, or just one, see Figure 5.
The computational strategy is the following: • For a fixed biofilm configuration, we let the diffusion problems for C o , C a and ε relax to stationary values, which occurs in a short time scale τ Each tile contains cells of certain volume and energy, with an associated EPS volume around it, whose dynamics is governed by the DEB equations. Here, we assume that each tile contains just one cell. In practice, all these items are stuck together in a more amorphous mass. This is a simple computational idealization allowing us to combine macroscopic and microscopic information. • We revise the status of all cells in the biofilm, in a random ordering. Cells differentiate and divide according to the rules stated in Section 2.4. To decide whether a cell dies, we generate a random number n ∈ (0, 1) and kill the cell if p < n Na N , where N is the total number of cells and N a the current number of alive cells.
• We update the concentration fields and repeat the procedure.
Allocating newborn cells, reabsorbing dead cells and allocating produced EPS in a 2D or 3D biofilm geometry are challenging issues for which many Each tile contains cells of certain volume and energy, with an associated EPS volume around it, whose dynamics is governed by the DEB equations. Here, we assume that each tile contains just one cell. In practice, all these items are stuck together in a more amorphous mass. This is a computational idealization allowing us to combine macroscopic and microscopic information.
ideas have been proposed depending on the biofilm types [10,11,12,24,25,39]. We do not intend to address these issues here. Instead, we work in a simplified geometry, implementing simple rules that allow us to check the performance of the DEB dynamics in a spatially varying film. In the numerical simulations presented here we make the following assumptions: • The biofilm is a 2D region, divided in cubic tiles for computational purposes.
• Each tile contains one cell together with EPS and dissolved substances. The idea is schematically represented in Figure 5.
• When a cell divides, the length and energy split in two. The daughter cell shifts other cells in the direction of minimal mechanical resistance (closest to the surface or to a dead cell).
• When a cell dies, it is reabsorbed by its neighbors, provided there are enough neighbors alive. A neighbor or newborn cell occupies its place in that case. Otherwise, its stays there and we get a necrotic region.
• A fraction α ∈ (0, 1) of the EPS produced by equation (3) remains attached to the cell. The rest diffuses through the whole biofilm. Therefore, when updating the concentrations C a , C o , ε, we need to update the diffused concentration of EPS C e (x, t) too: with zero flux boundary conditions ∂Ce ∂n = 0 on the biofilm boundary. The equations in Sections 2.2-2.3 set α = 1. All the EPS produced by a cell remains attached to it. When α = 1, R e is replaced by αR e + R e in all these equations. αR e represents the fraction of EPS produced at each site that remains in it. It may be zero when there are no EPS producers in that tile. R e represents the fraction of EPS accumulated due to global diffusion processes governed by (14). It will always be non zero by diffusion. The way we evaluate R e is the following. If at a certain control volume we start with an initial concentration C e (t 1 ) to reach a final concentration C e (t 2 ) this defines a rate R e = 1 t 2 −t 1 ln( Ce(t 2 ) Ce(t 1 ) ). The overall computational chart is represented in Figure 6.
The physical justification for this treatment of EPS is that the EPS matrix is formed by polymers. Monomers will diffuse easily, but a fraction of them will form polymeric chains of increasing size attached to the cells. An empirical justification is that if the produced EPS remains attached to producers, the active cells not producing EPS in outer layers die too fast and form unphysical necrotic layers in absence of any toxicants. In practice, the EPS matrix envelops and shelters all the biofilm cells.  is defined as the distance into the biofilm at which the first derivative of the oxygen concentration reaches 5% of its maximum value. As said before, oxygen is considered to be the limiting factor [16] for biofilm development in medical environments. We lack measurements of k, k , Y x/o and Y p/o in this case, but we may propose values for Y x/o and Y p/o , k and k based on an educated guess respecting the proportions observed for carbon [34].

Parameter calibration
Values for the maintenance rates m κ and investment ratios g κ are inferred as follows. The macroscopic equations for EPS and cell production are [34]: In absence of acclimation, toxicity and harzard corrections, that is a = a M , ν = ν, h = 0, equations (2) and (3) read Notice that the energy governed by (1) tends to an equilibrium e = f for each cell. We set f = f max = Co,out Co,out+Ko and e = f max . Dimensions can be restored in these equations by scaling v and v e . Comparing with the equations for X e and X, we approximate gκ mκ by k k and νfmax−mκgκ fmax+gκ by µ max f max . Taking m κ = k k g κ , g κ is a positive solution of a second order equation, given by The values of the multiplicative stress coefficient s q , Weibull aging acceleration h a , dissolved antibiotic toxicity k I q,A , influx and efflux coefficients k I A , k O A , and the initial hazard rate h 0 in equations (10), (11), (12) and (13) are fitted to data in reference [17]. To do so, we consider a single cell of volume v, energy e with density ρ x and write down the equations for the hazard rate h, aging acceleration q and survival probability p set- .
The way we find values for a, b, c, d, s is the following. We first remove the antibiotic and seek to adjust the viable cell counts in Table 1 of reference [17] for the control population. Figure 7 compares the bacterial counts at times 3, 6, 24 hours with our fittings of log 10 (pN ) for a population of the initial size N = 10 5.75 CFU/m. Once these parameters are fitted, we successfully fit the remaining parameters to the approximate number of viable cells for four different antibiotics: levofloxacin, ciprofloxacin, gentamicin, ceftazidime.
We approximate the data in Table 1 of [17] in absence of antibiotics dividing by 10 the values s q and h a from reference [21]. We choose average values for e and v, v = 0.5, e = median(f (C o,out )) ≈ 0.003, and set ν = ν. Then, c = e h a (ν − r) ≈ 3.61 × 10 −7 y d = e s q ρ x v(ν − r) ≈ 4.8 × 10 −3 . We vary a and b so as to approximate the data in Table 1 of [17]: a = h 0 = 0.4 and b = 0.29. These values are chosen because the fit still holds when we replace single cells by a biofilm formed by a few hundred layers. Since r ∼ 0 for EPS cells, we set k = b and keep the ratio k k observed for carbon.

Numerical results
The purpose of this section is to investigate the influence of the spatial variations of the concentrations in cell death. To simplify, in the small scale 2D simulations presented here we consider tiles of size 1 µm × 1 µm in Figure  5. Figure 8 shows the spatial structure of the oxygen distribution in a biofilm of maximum width 1000 µm and maximum height 200 µm. We solve equation (4) discretized in the grid, with Dirichlet boundary condition C o = C o,out on the interface with air and Neumann boundary condition dCo dn = 0 on the bottom interface. We take f (C o ) = Co Co+Ko and X = ρ x v, where v is the dimensionless volume of the cell contained in the tile. Figure 8 sets v = 0.5. Initially, C o = 0 inside the biofilm. We solve using a explicit scheme until C o Figure 7: Comparison of theoretical predictions and experimental data. The predictions provided by formula (15) for p are compared to the data in Table 1 of Reference [17]. The parameter values for p(t) are a = 0.4, b = 0.29, c = 3.61 × 10 −7 and d = 4.8 × 10 −3 . The additional parameter s is adjusted for each antibiotic. Circles represent log 10 (p · N ) with N = 10 5.75 (magenta) and N = 10 5.75±0.10 (red). Asterisks represent the measurements in Table 1 of Reference [17] for the same antibiotic doses. As before, light (blue) colors correspond to the mean measured value, whereas dark (black) colors mark the error interval.  relaxes to a stationary configuration: where g(X) represents the source in equation (4). The spatial step δx is the side of a tile. Tiles in the spatial grid are labelled using a couple of indices (i, j) to indicate their position. Time is discretized with time step δt, so that t +1 = t + δt = ( + 1)δt and C o;i,j ∼ C o;i,j (t ). For this scheme to be stable, we need δt do δx 2 ≤ 0.5. We advance in a time scale of seconds. Similar schemes are used to update the additional space dependent concentrations C a , C ε , C e when needed.
The cell energies are governed by DEB equations depending on f = Co Co+Ko , which varies through the biofilm. Figure 9 represents the equilibrium distribution e = f (C o ) at each tile. The maximum energy is e m = Co,out Co,out+Ko . In view of the energy distribution, we set threshold values e eps1 = 7 × 10 −6 , e eps2 = 0.02 and e d = 0.01 to produce EPS or divide.
Starting from a uniform distribution of cell types (non EPS producers, v e = 0) with volumes v(0) = 0.5 we compute the evolution of their energy and size to check the evolution of the cell type distribution with time. We choose as initial energy a perturbation of the equilibrium energy distribution at each location, and fix q > q eps = 10 −8 . By solving the equations for e, v, v e for each cell tile ignoring aging and hazard, we see that cells evolve to be EPS producers in an intermediate region. The outer layers are normal active cells, likely to divide. A thin inner core may be formed by cells with little energy than do not divide and do not produce EPS. The equilibrium distribution of types is depicted in Figure 10(a). If we increase the initial energy, we may find intermediate states such those in Figure 10(b), where that thin inner layer vanishes. Repeating the simulation with initial volumes v randomly distributed between (0.45, 0.6) and random perturbations of the equilibrium energy we find similar results with a smoother transition between regions, reflecting the noise.
To study the distribution and number of dead cells, as well as the effect of antibiotics, we activate the equations for aging and hazard rates. We compute a reference concentration of antibiotic C a in the biofilm solving (5) with Dirichlet boundary condition C a = C a,out on the interface biofilm/air and Neumann boundary condition dCa dn = 0 on the bottom interface. We assume that all cells are undifferentiated in this reference computation (R e = 0) and employ the equilibrium oxygen and energy distributions to evaluate R. The concentration of antibiotics inside the biofilm evolved towards a constant, and we will take it to be approximately constant and equal to C a,out in the biofilm. Figure 11 shows that the fittings of the time evolution of the survival rates performed in Section 3 persist when all the biofilm is taken into account. A low energy cell population is considered in that figure, forbidding cell division for simplification.
The model yields information on qualitative behavior that agrees with Figure 11: Simulated live cell counts compared to experimental measurements in Table 1 of Reference [17]: (a) in absence of antibiotics, (b) for ceftazidime. The continuous curve represents the simulated values. Asterisks depict the measurements in Table 1 of reference [17] for the same antibiotic doses, black lines represent the error bars.
some experimentally observed tendencies and may be useful to understand the causes or to predict new behaviors. Panels (a) and (c) in Figure 12 reveal that for an antibiotic acting mostly on active cells like Ceftazidime (it inhibits cell wall synthesis), cells mostly die in the outer biofilm layer and the extent of the dead region increases with the dose of antibiotic. Panel (b) shows that it also increases with time. We are assuming that the antibiotic has diffused inside the biofilm so that the spatial distribution of dead regions depends on what type of cells are preferential targets of the antibiotic under study. For antibiotics like Ceftazidime, acting on active cells, we use toxicity coefficients depending on the cell activity, of the form e 0.1 median(e) k I qA . For other antibiotic like colistin, targeting the dormant core cells, we may use toxicity coefficients decreasing with the cell activity instead. Mutations increasing the efflux might be accounted for raising the value of that coefficient. Augmented resistance to deactivating enzymes may be incorporated modifying the toxicity.
In absence of antibiotics cells die in the inner regions in a scattered way. Active cells in the outer layers are still more likely to die, a standard fact in the DEB theory, that assumes organisms with more energy more likely to suffer damage. That effect would be diminished increasing the difussivity of produced EPS d e , or allowing differentiation into EPS producers for larger energies. Cells that die in the outer layers in absence of antibiotics may also be eroded by external agents or absorbed by newly born cells. Unlike the previous case, we do not get an expanding outer necrotic region. Analyzing the amount of EPS produced with and without antibiotics, we find that the presence of antibiotics enhances EPS production, see Figure 13. This fact is also responsible for the small difference between panels (a) and (c) of Figure  12. The number of dead cells is larger for panel (c), but not as large as could be expected from the dose increase due to enhanced EPS production that neutralizes that fact.

Conclusions
Coupling dynamic energy budget dynamics for the metabolism of individual cells to diffusion equations for the evolution of chemical fields we have reproduced qualitative trends experimentally observed in bacterial biofilms under the action of antibiotics. We calibrate the parameters for a quantitative agreement with the death rates measured for P. aureginosa under a number of antibiotics and study spatial variations in the distribution of (a) (b) dead and alive cells. Our numerical simulations show that the presence of antibiotics enhances production of polymeric matrix, as expected from true measurements. For antibiotics targeting active cells, we observe the formation of a necrotic outer layer, that expands deeper in the biofilm as time goes on and the antibiotic dose increases. Our tests suggest ways to handle antibiotics targeting dormant cells and resistance mechanisms such as efflux pumps or enzymes. New experimental measurements [31] would be required to calibrate the mechanisms to practical cases. The current computational model enables validation and hypothesis testing when developing therapies to handle chronic infections caused by biofilms. Determining the true relevance of increasing the dose and whether periodic or continuous infusion of antibiotics is more effective [20] would be essential for practical applications.


Single-cell transcriptomics identifies CD44 as a new marker and regulator of haematopoietic stem cells development

The endothelial to haematopoietic transition (EHT) is the process whereby haemogenic endothelium differentiates into haematopoietic stem and progenitor cells (HSPCs). The intermediary steps of this process are unclear, in particular the identity of endothelial cells that give rise to HSPCs is unknown. Using single-cell transcriptome analysis and antibody screening we identified CD44 as a new marker of EHT enabling us to isolate robustly the different stages of EHT in the aorta gonad mesonephros (AGM) region. This allowed us to provide a very detailed phenotypical and transcriptional profile for haemogenic endothelial cells, characterising them with high expression of genes related to Notch signalling, TGFbeta/BMP antagonists (Smad6, Smad7 and Bmper) and a downregulation of genes related to glycolysis and the TCA cycle. Moreover, we demonstrated that by inhibiting the interaction between CD44 and its ligand hyaluronan we could block EHT, identifying a new regulator of HSPC development.

Introduction
Understanding the developmental origin of haematopoietic stem and progenitor cells (HSPCs) is of critical importance to efforts to produce blood and blood products in vitro for medical applications. Haematopoietic stem and progenitor cells (HSPCs) originate from endothelial cells in the aorta gonad mesonephros (AGM) of midgestation embryos 1,2 . This process known as the endothelial to haematopoietic transition (EHT) requires drastic morphological changes that have been directly visualised through time-lapse imaging studies both in vitro and in vivo [3][4][5][6] . EHT is a highly conserved process that has been studied across vertebrate models from xenopus and zebrafish to mice 7 . Importantly, the human definitive blood system also has an endothelial origin 8 .
The best tools so far to detect endothelial cells with haemogenic capabilities rely on using fluorescent reporters under the control of Runx1 9 or Gfi1 10 regulatory elements, two keys transcription factors in the process of EHT. Cells expressing these transcription factors already co-express blood and endothelial genes. However, we still do not know the nature of endothelial cells which will acquire the expression of The early haematopoietic hierarchy has been described as a three-step process based on phenotypic characteristics. Specifically, Pro-HSC, Pre-HSC type I and type II populations have been defined based on their expression of the cell surface markers CD41, CD43 and CD45, as well as the time taken mature into definitive HSCs in OP9 co-culture 13 . Recently, an in depth transcriptional investigation was performed on the 4 type I and type II Pre-HSC populations at day 11 of mouse embryonic development, however, the earlier stages of EHT remain largely uncharacterised 14 . Indeed, gaining a solid understanding of HE and the initial steps that endothelial cells must take to become HSPCs has proved difficult in the absence of a robust marker.
Through antibody screening and single-cell RNA sequencing (sc-RNA-seq) we identified CD44 as a novel marker of EHT, enabling the isolation of key cellular stages of blood cell formation in the embryonic vasculature. CD44 is a cell surface receptor principally involved in the binding of the extracellular matrix molecule hyaluronan 15 . Its cell surface expression has been used to identify cancer stem cell populations and has been strongly linked to the metastatic potential of many cancers [16][17][18][19][20] . While previous research has revealed the importance of CD44 in HSPCs migration to the bone marrow, the role of the receptor in early embryonic haematopoiesis has not been characterised 21 . Using CD44 expression in conjunction with VE-Cadherin (VE-Cad) and Kit we could clearly differentiate between vascular endothelium, HE, Pre-HSPC-I and Pre-HSPC-II more accurately than using the combination of VE-Cad, CD41, CD43 and CD45 markers. This has allowed us to perform extensive transcriptional profiling making it possible to characterise the very earliest changes in haematopoietic differentiation from endothelial cells. Moreover, by disrupting the interaction of CD44 and its ligand, we could inhibit EHT, demonstrating an unexpected role for CD44 in the emergence of HSPCs.

Antibody screening and single-cell RNAseq identifies CD44 as a potential marker of early haematopoietic fate
It is well established that HSPCs have an endothelial origin within the embryo [4][5][6]22 . In order to better characterise the transition between endothelial and haematopoietic identity, we performed both an in vitro antibody screen and an in vivo single-cell RNAseq experiment to identify new markers allowing the identification of subpopulations within VE-Cad + cells ( Supplementary Fig. S1). Antibodies against 176 cell surface markers 23 were tested against the VE-Cad + population generated from our in vitro ESC differentiation system into blood cells which recapitulates faithfully the early blood development process 24 . Forty-two of these markers were expressed on VE-Cad + cells (Supplementary Table S1). We looked for bimodal expression to separate distinct endothelial populations and identified a short-list of sixteen candidates to test in vivo including CD41 and CD117 (also known as Kit, a marker of HSPCs) already known to split VE-Cad + cells 13 . Eight of these markers (CD44, CD51, CD55, CD61, CD93, Icam1, Madcam1 and Sca1) were found to split the VE-Cad + endothelial population of the AGM in two (Fig. 1a). In parallel, VE-Cad + cells were isolated from the AGM region of E10.5 embryos and their transcriptional profiles analysed sc-RNA-seq ( Fig. 1b-d). Clustering analysis identified a population with both haematopoietic and endothelial gene expression, distinct from the other endothelial population (Fig. 1d). Bioinformatics analysis showed that Cd44 is one of the best marker genes for this population of transitioning cells co-expressing endothelial and haematopoietic genes (Fig. 1c). The expression of Cd44 was also positively correlated with other known haematopoietic markers such as Runx1, Gfi1 and Adgrg1 (Gpr56) (Fig. 1d).
Given the association of Cd44 with endothelial cells undergoing EHT at both the protein and mRNA level we decided to further investigate its role in embryonic haematopoiesis.

CD44 marks transitioning cells with differing morphology
To validate our screening results and investigate the identity of CD44 + cells, we performed immunofluorescence and more detailed flow cytometry analysis on the 6 AGM region of mouse embryos at embryonic day 9.5 (E9.5) and 10.5 (E10.5) (Fig. 2).
Immunofluorescence of cross-sections of mouse AGMs revealed that CD44 marked cells that were part of the vascular wall and cells that were incorporated in haematopoietic clusters (Fig 2a). Flow cytometry revealed that CD44 expression significantly increased in the VE-cad + endothelium of the AGM between E9.5 and E10.5 when cells are undergoing EHT (Fig. 2b-c). Furthermore, by staining with an antibody against Kit (a marker of intra-aortic haematopoietic clusters), we found that the majority of cells with lower levels of CD44 expressed little or no Kit (Fig 2d).
By grouping the cells based on their cell surface expression of CD44 and Kit, we found these populations to be significantly different in terms of cell size, suggestive of cells undergoing a morphological transition (Fig. 2e). Altogether, our results indicate that CD44 marks a subset of endothelial cells and cells in the haematopoietic clusters of the AGM during the key window of HSPCs development in the mouse embryo.

Single-cell q-RT-PCR analysis identifies three homogeneous CD44 + populations with an increasing haematopoietic profile
Using the Biomark HD single-cell qPCR platform, we analysed the expression of 95 genes associated with both endothelial and haematopoietic cell types 24 . We performed extensive transcriptional profiling on the CD44 Neg , CD44 Low Kit Neg , CD44 Low Kit Pos and CD44 High populations identified (Fig. 2d) between E9.5 and E11.5 (342 cells in total). We found that the CD44 Low Kit Pos and CD44 High populations expressed numerous haematopoietic genes (Fig. 3a). The CD44 Low Kit Pos cells also expressed many endothelial genes as in our sc-RNA-seq analysis (Fig. 1d), however the CD44 High cells appeared to be more advanced in the EHT process and lacked endothelial gene expression (Fig. 3a). Of note, the CD44 Low Kit Pos population expressed specifically Gfi1 and Itgb3 (coding for CD61) as well as the heptad of transcription factors (Gata2, Runx1, Lyl1, Erg, Fli1, Lmo2 and Tal1) whose simultaneous co-expression is responsible for the dual endothelial-haematopoietic identity of Pre-HSPCs ( Supplementary Fig. S2) 24 .
Conversely, the CD44 Neg and CD44 Low Kit Neg populations both showed specific endothelial gene signatures and lacked haematopoietic gene expression (Fig. 3a).
Despite their different CD44 expression patterns, these cell populations clustered together ( Supplementary Fig. S2). We repeated this experiment using a new selection 7 of 96 genes based on our sc-RNA-seq experiment (Supplementary Table S2). With our new gene list, we were able to confirm the dual endothelial-haematopoietic and haematopoietic identities of the CD44 Low Kit Pos and CD44 High populations, respectively. Surprisingly, the CD44 Neg and CD44 Low Kit Neg populations formed two distinct clusters (Fig. 3b). In addition to Cd44, we found the genes Adgrg6, Emb, Fbn1, Pde3a, Plcg2, Serpinf1, Smad6, Smad7, Sox6 and Stxbp2 to be up-regulated in the CD44 Low Kit Neg cells compared to CD44 Neg . In contrast, Bmp4, Kit, Hmmr and Pde2a were more expressed in CD44 Neg endothelial cells (Fig. 3c).
Although the four groups defined by VE-Cad, CD44 and Kit were confirmed to be distinct transcriptionally, our clustering analysis found a fifth population (SC_3) composed of cells from both CD44 Low Kit Neg and CD44 Low Kit Pos groups (Fig. 3a). Its transcriptional profile was found to be intermediary between SC_2 (CD44 Low Kit Neg ) and SC_4 (CD44 Low Kit Pos ) e.g. it expressed Adgrg1, Runx1, Itgb3 and Spi1 like SC_4 but was still expressing Adgrg6 and Pcdh12 like SC_2 (Fig. 3c).
This suggests a developmental link between the CD44 Low populations where CD44 Low Kit Neg cells would be the direct precursors of the CD44 Low Kit Pos population which would then go on to generate CD44 High cells. Interestingly, it is within this transitional SC_3 population that we saw the up regulation of Runx1, Spi1 and Gfi1, which are three of the four transcription factors used to reprogram adult mouse endothelial cells into HSCs 12 .
VE-Cad, CD44 and Kit could segregate the earliest stages of haematopoietic development more accurately than the combination of VE-cad, CD41, CD43 and

CD45 markers
To place our results in context of known AGM subpopulations, we performed transcriptional analysis of the Pro-HSC, Pre-HSC type I, and Pre-HSC type II populations defined according to the combination of VE-cad, CD41, CD43 and CD45 markers (Fig. 4a) 13 . Following hierarchical clustering, we found three clusters: the first mostly composed of Pro-HSCs, a second being a mix of Pro-HSCs and Pre-HSCs type I and a third one composed only of Pre-HSCs type II ( Supplementary Fig. S3).
We then performed a clustering analysis in conjunction with the populations defined by CD44 (Fig. 4b). This revealed that the SC_5 (CD44 High ) population closely associated with the Pre-HSC type II and the SC_4 (CD44 Low Kit Pos ) population with Pre-HSC type I and part of the Pro-HSCs. SC_2 (CD44 Low Kit Neg ) clustered closely 8 with the remaining Pro-HSCs. Finally, only five cells with Pre-HSCs type I and Pro-HSCs phenotype clustered with the SC_1 (CD44 Neg ). Ninety-seven per cent of Pro-HSCs, Pre-HSCs type I, and Pre-HSCs type II were CD44 positive ( Supplementary   Fig. S3 and Fig. 4c).
Overall, we have demonstrated that the phenotypes based on CD44 expression could allow us to isolate all key populations in the process of HSCs formation more accurately than the phenotypes previously described.
Bulk RNA-seq analysis further distinguishes CD44 Neg and CD44 Low Kit Neg endothelial populations and identifies early changes in the differentiation process To further compare the two endothelial populations found in the AGM, we performed RNAseq on 25-cell bulk samples from these populations across three different time points (E9.5, E10 and E11). We analysed as well the more advanced stages in EHT: CD44 Low Kit Pos at E9.5 and E10 and CD44 High at E11. The bulk RNAseq approach allowed us to detect low abundant genes (such as genes encoding transcription factors) more efficiently than sc-RNA-seq and also to measure smaller changes in gene expression between populations.
The samples clustered according to their marker expression, despite the difference in developmental time, confirming our previous experiments (Fig. 5a). We identified several haematopoietic genes switched on in the CD44 Low Kit Neg population including, Ctsc, Nfe2, Runx1 and Ifitm1, suggesting that these cells are subjected to the EHT process (Fig. 5b). The expression pattern of endothelial genes fits with our previous observations -more highly expressed in CD44 Neg and CD44 Low Kit Neg , moderately expressed in CD44 Low Kit Pos and absent in CD44 High (Fig 5b). Moreover, we used this dataset to check the expression of genes corresponding to proteins we found in our antibody screen (Fig. 1a). Cd93 and Madcam1 showed an endothelial expression pattern and could be used to separate endothelial from blood cells in the AGM. In contrast, Itgb3 marks specifically the CD44 Low Kit Pos population while Ly6e marks both CD44 Low Kit Pos and CD44 High populations as shown previously in Supplementary Fig. S2 and Fig. 3.
Interestingly, this transcriptome analysis showed strong differences between the two endothelial populations of the AGM. We found 1605 genes differentially expressed between these two populations (p-value < 0.01, Wald test). Among them, several genes from the Notch pathway (Hey2, Jag1, Dll4, Hey1 and Notch1) were 9 significantly more expressed in the CD44 Low Kit Neg population compared to CD44 Neg .
Similarly, antagonists of the TGFbeta/BMP pathway, including Smad6, Smad7 and Bmper, were up-regulated in CD44 Low Kit Neg cells compared to CD44 Neg ones. In contrast, target genes of the Wnt pathway (Lef1, Ccnd1 and Myc) were more highly expressed in CD44 Neg compared to CD44 Low Kit Neg . From this analysis we could also identify specific markers for each of the endothelial populations (Fig. 5c). Remarkably, a large number of the differentially expressed genes between the two endothelial populations belonged to metabolic processes (395 out of 1605 genes; 1.32-fold enrichment, p-value <0.05, Fisher's exact test). We therefore identified specific metabolic pathways distinguishing the two populations ( Fig. 6a and Supplementary Table S3). Notably, the CD44 Low Kit Neg population showed a pronounced down regulation of genes coding for enzymes involved in glycolysis, TCA cycle and respiration, suggesting reduced ATP generation. Furthermore, amino acid and nucleotide biosynthesis genes were also down regulated. Altogether it suggests that the CD44 Low Kit Neg is a non-proliferative, metabolically rather quiescent state which is in line with the smaller size of this population compared to the CD44 Neg 10 ( Fig. 2e). Given that endothelial cells are known to obtain most of their energy from glycolysis 25 , this change in metabolic status supports a loss of endothelial identity.
These cells also show a marked increase in the expression of pathways leading to lipids with regulatory function: glycerolipids, glycerophospholipids, phosphatidylinositol and sphingolipids (Fig. 6a).
Interestingly, we observe that several genes involved in autophagy, a process known to be regulated by phosphatidylinositol and sphingolipids 26-29 , are also highly upregulated in the CD44 Low Kit Neg population (Fig. 6b). Concordantly, two key processes accompanying autophagy, ubiquitylation and proteolysis, are also upregulated. As autophagy has been shown to play a key role in embryonic development and haematopoiesis [30][31][32] , this provides further support to the CD44 Low Kit Neg cells being in transit from endothelial to haematopoietic cells.

Runx1 is not required for the formation of CD44 Low Kit Neg endothelial cells
CD44 has allowed us for the first time to clearly define the key VE-Cad + populations in the AGM. The transcription factor RUNX1 is a key driver of HSPC development and is known to down-regulate endothelial identity through its target genes GFI1 and GFI1B 10,33 . Next, we decided to evaluate the impact of Runx1 loss of function on the different CD44 + cells. Using a Runx1 knockout mouse model 34 , we stained for VE-Cad, CD44 and Kit expression and performed transcriptional profiling on the sorted cells (Fig. 7). We found that in the absence of Runx1 there is a loss of CD44 High and CD44 Low Kit Pos cells and we observed a concomitant increase in the frequency of the CD44 Low Kit Neg population ( Fig. 7a and 7b). Interestingly, we found no obvious transcriptional differences between the CD44 Low Kit Neg populations derived from Runx1 +/+ versus Runx1 -/embryos indicating that Runx1 is not necessary for the formation of these endothelial cells but for the promotion of the transition into CD44 Low Kit High cells ( Supplementary Fig. S5).

All CD44 + populations display haematopoietic potential
To understand the haematopoietic potential of the different populations defined by VE-Cad, CD44 and Kit expression, we performed ex vivo assays using an OP9 co-culture system. No colonies were formed at the single-cell level from either the CD44 Low Kit Neg or the CD44 Neg populations. However, by plating cells at a density 11 of 300 cells per well, we could observe round cell colony formation from the CD44 Low Kit Neg population but not from the CD44 Neg (Fig. 8a).
In contrast, using single-cell sorting, we found that the CD44 Low Kit Pos population was the most potent with an average of 43% of single cells forming round cell colonies after three days of growth. Similarly, the CD44 High population produced haematopoietic colonies but with a lesser frequency, on average 11% of single cells showed the ability to form haematopoietic colony on OP9 (Fig. 8b-c). To uncover the differentiation potential of the cells generated by CD44 Low Kit Pos and CD44 High , we performed CFU assays following the OP9 co-culture. Both populations readily generated both erythroid and myeloid colonies with the CD44 Low Kit Pos population demonstrating a significantly higher capacity than the CD44 High population ( Fig. 8d-e).
However, the four-fold increase in the number of CD44 Low Kit Pos colonies on OP9 did not correspond to a four-fold increase in CFU colonies suggesting a higher replating efficiency of the CD44 High cells compared to CD44 Low Kit Pos (Fig 8f). We further tested the lymphoid potential of the CD44 High population by growing cells for 21 days on either OP9 or OP9-DL1 with lymphoid promoting cytokines. We demonstrated that this population could give rise to both B and T cells ex vivo (Fig. 8h).
Overall, we found that all populations expressing CD44 displayed haematopoietic differentiation capacity including CD44 Low Kit Neg reinforcing the differentiation link between them as suggested by the transcriptome analyses described previously.

Interrupting hyaluronan binding to the CD44 receptor inhibits the endothelial to haematopoietic transition ex vivo and in vitro
So far, we demonstrated that CD44 was a very useful marker to distinguish the different stages of EHT. Although the CD44 knock-out mice do not have a severe haematopoietic phenotype 35 , compensatory mechanisms through other Hyaluronan receptors (e.g. Hmmr 36 expressed by some CD44 Low Kit Pos and CD44 High cells in Fig.   3a) may be at play to diminish the consequences of CD44 loss of function. In order to explore the functional role of CD44 in EHT we employed a pharmacological approach. By treating CD44 High sorted cells with a CD44 blocking antibody, we found that round cell colony formation could be inhibited in a dose dependent manner ( Fig.   9a-b). The blocking antibody inhibited not only the number of colonies deriving from the ex vivo sorted cells but also the size of the colonies generated (Fig 9c). 12 To further investigate the role of CD44 in EHT we used an ESC in vitro differentiation system that mimics embryonic haematopoiesis. We performed a haemangioblast culture and analysed CD44 expression at day 1, 2 and 3. Endothelial In conclusion, these results for the first time demonstrate a regulatory role for hyaluronan and its receptor CD44 in the formation of HSPCs. 13

Discussion
Using antibody screening and sc-RNA-seq to dissect the endothelial populations in the AGM, we discovered that CD44 was a robust marker to distinguish all the main populations in the EHT process. Combining CD44 with Kit and VE-Cad allowed us to discriminate the different stages of EHT more accurately that the method based on VE-Cad, CD41, CD43 and CD45 cell surface markers 13 . In addition, we showed that CD44 had an unexpected regulatory function in the EHT process.
Our work has been instrumental in distinguishing the different types of endothelial cells in the AGM region. The CD44 Low Kit Neg has a gene expression signature strongly compatible with arterial identity (e.g. expression of Efbn2 and Sox17 and upregulation of Notch pathway target Hey2) while the CD44 Neg cells coexpressed genes related to the venous (e.g. Nr2f2 and Aplnr) and arterial cell fates (e.g. Efnb2 and Sox17) (Fig. 5d). This co-expression of venous and arterial genes supports previous work indicating that the dorsal aorta can contribute to the formation of the cardinal vein 37 .
Another interesting finding was the striking metabolic state difference between the two endothelial populations. The CD44 Neg cells appeared much more metabolically active than the CD44 Low Kit Neg suggesting that the latter population is becoming quiescent. It is surprising since the acquisition of a quiescent phenotype in endothelial cells occurs normally after birth following the completion of angiogenesis. Smad7 at the single-cell level suggesting that the quiescence phenotype we observed could also be linked to the co-expression of the two inhibitory Smads.
The acquisition of quiescence in this context could be the first indication of the divergence from an endothelial phenotype. The metabolic distinction of the CD44 Low Kit Neg population was coupled with an increase in genes involved in 14 autophagy. This is in line with the known role of autophagy in embryogenesis, haematopoiesis and stem cell maintenance 30-32,40,41 , such as protein and organelle turnover and protection from reactive oxygen species. Our data thus suggest that these cells mark the resetting of metabolic and regulatory states to initiate the EHT process.
We consequently propose that the CD44 Low Kit Neg population is the source of HSPCs in the AGM, hence the haemogenic endothelium ( . Given the role of CD44 in the metastatic process, it is possible that there is an overlap in its function for these transformations. Therefore, understanding the downstream targets of the CD44-hyaluronan interaction could also have implications also for cancer biology. 15

Timed mating and embryo dissection
For timed pregnancies, C57BL/6 wild-type mice or Runx1 +/mice were mated overnight. Embryos were collected in PBS supplemented with 10% FBS (PAA Laboratories

Single-cell RNA sequencing data analysis
The paired 2 x 101bp Illumina reads from the libraries were quantified using Salmon 43 with the setting -l IU to indicate library topology, and the optional flags --posBias and --gcBias to account for coverage and amplification biases present in sc-RNA-seq protocols. As an index the cDNA annotation of Ensembl release 85 for GRCm38.p4 was used, together with ERCC spike-in sequences. The TPM values were rescaled to not include ERCC expression and only consider endogenous gene expression.
Technical features of the data were compared with manual annotation of samples in C1 chambers through microscopy. Samples with less than 500,000 mapped reads and more than 30% mitochondrial content were discarded from analysis. This Novel markers for hematopoietic cells were identified using a likelihood ratio test, where the alternative model included a binary term for whether the cell was haematopoietic, and the null model just assumed a common mean for all the cells.
The P-values from the likelihood ratio test were corrected for multiple testing by the Bonferroni procedure. The top differentially expressed genes were investigated to find markers which could be used in follow-up experiments, and Cd44 was considered a good candidate (Fig. 1c).

In vitro ES cell differentiation system
The

Antibody screen, Flow cytometry and cell sorting
The antibody screen was performed using the Mouse Cell Surface Marker  imaged on a Leica SP5 confocal microscope.

Single-cell qPCR data analysis
Analysis of single-cell qPCR data was performed as previously described 24 .
Briefly, initial analysis was performed using the (Abcam). The number of colonies were assessed after 3 days in culture.

Haematopoietic colony forming assay
One hundred cells were initially sorted onto a confluent OP9 stromal layer as per OP9 co-culturing assay. After three days in culture cells were harvested with TrypLE express (Gibco) and colony-forming unit-culture (CFU-C) assays were initiated using Methocult complete medium (Stem Cell Technologies). Cells were grown in 35mm culture dishes and colonies quantified after 7 days.  M  o  r  g  a  n  O  a  t  l  e  y  ,  C  o  n  c  e  p  t  u  a  l  i  z  a  t  i  o  n  ,  F  o  r  m  a  l  a  n  a  l  y  s  i  s  ,  I  n  v  e  s  t  i  g  a  t  i  o  n  ,  V  i  s  u  a  l  i  z  a  t  i Figure S1 and File S1.        sc1_cell1_salmon_out  sc1_cell71_salmon_out  sc1_cell53_salmon_out  sc1_cell61_salmon_out  sc1_cell84_salmon_out  sc1_cell90_salmon_out  sc1_cell10_salmon_out  sc1_cell58_salmon_out  sc1_cell92_salmon_out  sc1_cell17_salmon_out  sc1_cell2_salmon_out  sc1_cell25_salmon_out  sc1_cell36_salmon_out  sc1_cell34_salmon_out  sc1_cell12_salmon_out  sc1_cell65_salmon_out  sc1_cell32_salmon_out  sc1_cell52_salmon_out  sc1_cell87_salmon_out  sc1_cell86_salmon_out  sc1_cell20_salmon_out  sc1_cell57_salmon_out  sc1_cell38_salmon_out  sc1_cell68_salmon_out  sc1_cell50_salmon_out  sc1_cell56_salmon_out  sc1_cell22_salmon_out  sc1_cell33_salmon_out  sc1_cell55_salmon_out  sc1_cell51_salmon_out  sc1_cell19_salmon_out  sc1_cell4_salmon_out  sc1_cell63_salmon_out  sc1_cell94_salmon_out  sc1_cell47_salmon_out  sc1_cell48_salmon_out  sc1_cell96_salmon_out  sc1_cell23_salmon_out  sc1_cell78_salmon_out  sc1_cell79_salmon_out  sc1_cell27_salmon_out  sc1_cell67_salmon_out  sc1_cell44_salmon_out  sc1_cell5_salmon_out  sc1_cell62_salmon_out  sc1_cell80_salmon_out  sc1_cell43_salmon_out  sc1_cell69_salmon_out  sc1_cell64_salmon_out  sc1_cell66_salmon_out  sc1_cell13_salmon_out  sc1_cell49_salmon_out  sc1_cell54_salmon_out  sc1_cell60_salmon_out  sc1_cell40_salmon_out  sc1_cell28_salmon_out  sc1_cell42_salmon_out  sc1_cell76_salmon_out  sc1_cell26_salmon_out  sc1_cell29_salmon_out  sc1_cell18_salmon_out  sc1_cell82_salmon_out  sc1_cell93_salmon_out  sc1_cell70_salmon_out  sc1_cell46_salmon_out  sc1_cell74_salmon_out  sc1_cell16_salmon_out  sc1_cell14_salmon_out  sc1_cell45_salmon_out  sc1_cell35_salmon_out  sc1_cell77_salmon_out  sc1_cell8_salmon_out  sc1_cell11_salmon_out  sc1_cell15_salmon_out  sc1_cell6_salmon_out  sc1_cell30_salmon_out  sc1_cell39_salmon_out


s2fieldsofstudy: Biology Medicine

text:
Association between AGTR1 A1166C polymorphism and the susceptibility to diabetic nephropathy

Abstract Background: Diabetic nephropathy (DN) is a common complication in patients with diabetic mellitus (DM). Growing evidences have demonstrated that the polymorphisms of angiotensin II receptor type 1 (AGTR1) showed significant association with DN onset, but no consensus has been achieved yet. Therefore, we performed this meta-analysis to combine the findings of previous researches for a more comprehensive conclusion. Methods: Eligible publications were identified through electronic databases. The intensity of the correlation between AGTR1 A1166C polymorphism and DN susceptibility was evaluated through calculating pooled odds ratios (ORs) with their 95% confidence intervals (95% CIs). Heterogeneity among included studies was examined with Q test. When P-value less than .05, significant heterogeneity presented, random-effects model was used to calculate the pooled ORs, otherwise, the fixed-effects model was used. Stratification analyses were also performed based on ethnicity and the type of DM. Results: Seventeen eligible articles were finally included in the present meta-analysis. The analysis results showed that AGTR1 A1166C polymorphism was significantly related to increased risk of DN under CC versus AA (OR = 1.723, 95% CI = 1.123–2.644), CC + AC versus AA (OR = 1.179, 95% CI = 1.004–1.383), CC versus AA + AC (OR = 1.662, 95% CI = 1.112–2.486), and C versus A (OR = 1.208, 95% CI = 1.044–1.397) genetic models. Additionally, a similar result was also found in Asian and T2DM (type 2 diabetic mellitus) groups after subgroup analyses of ethnicity and DM type. Conclusion: AGTR1 A1166C polymorphism may increase the susceptibility to DN, especially in Asians and T2DM population.

Introduction
Diabetic nephropathy (DN) refers to diabetic glomerulosclerosis, representing one of the frequently observed diabetes mellitus (DM) systemic microvascular complications. [1,2] According to statistics, the incidence rate of DN shows an increasing tendency in the past few years. [3] It is predicted by WHO that DM will be prevalent in developing countries in the 21st century. [4] Once persistent proteinuria occurs, it will progressively develop to endstage renal disease. [5] DN is one of the leading causes of disability and death in DM patients, being a critical topic in medical research. [6] Until now, the pathogenesis of DN could not be completely explained. Research data have shown that the occurrence of DN is related to multiple factors, such as changes in hemodynamics, metabolic disorders, and the involvement of growth factors and genetic elements. [7,8] While in recent modern medical studies, hereditary factors have been demonstrated to occupy an extremely vital position in the occurrence of DN. [9] For instance, a meta-analysis based on 1894 DN cases and 1746 controls demonstrated that NADPH oxidase p22phox C242T SNP showed obvious association with macroalbuminuria in patients with diabetes. [10] To investigate the genetic factors may provide a new insight into the pathogenesis of DN.
Renal hemodynamic abnormalities play an important role in the initiation and progression of DN. The abnormalities in reninangiotensin system (RAS), especially the local RAS of kidney, are the major cause of renal hemodynamic abnormalities. [11] Given the function roles of RAS in onset of DN, the alterations in RAS system genes might be involved in development of DN. [12] RAS is consisted of renin, angiotensin I-converting enzyme (ACE), angiotensinogen (AGT) and angiotensin II receptor, type 1 (AGTR1). [13] AGTR1 is widely expressed in various tissues, such as vessel walls, lung and kidney, and after being activated. AGTR1 can not only lead to water-sodium retention and elevated blood pressure but also participate in the microvascular disorders in type 2 DM (T2DM). [14] Meanwhile, activation of AGTR1 may regulate renal function. [15,16] The expression pattern of AGTR1 shows significant association with nephropathy. [17] There are several polymorphisms in the AGTR1 gene, including A1166C, T573C, A1062G, G1517T, and A1878G. Among them, the A1166C polymorphism is located at the 3 0 untranslated region of the gene, which does not affect the encoding process of AGTR1 protein in theory. But it still has the potential to influence the stability of the mRNA expression of the gene. [18,19] Growing evidences have proved the significant association between AGTR1 A1166C SNP and DN. [20][21][22][23] However, due to the differences in study ethnic population, type of DM, as well as the sample size, no conclusive result has been achieved yet.
In this study, we aimed to obtain a reliable result about the genetic association of AGTR1 A1166C polymorphism and DN susceptibility through a meta-analysis.

Study design
The present meta-analysis was performed based on the guidelines of the Preferred Reporting Items for Systematic Reviews and Meta-analyses (PRISMA) statement. The PRISMA checklist presented in the form of "Supplement information."

Literature searching
A systemic search was performed in the databases of PubMed, EMBASE, Google Scholar Web, China National Knowledge Infrastructure (CNKI), and Wanfang for eligible articles published in English or Chinese language, using the combination of the following key terms: "angiotensin II receptor, type 1" or "AT 1 receptor" or "AGTR1" or "AT1," "polymorphism" or "mutation" or "variant," and "diabetic nephropathy" or "DN" or "nephropathy." Besides, the reference lists of relevant articles were also manually checked for additional publications.

Selection criteria
All eligible articles had to satisfy the following criteria: with a case-control design; the individuals in control group were DM without DN patients, while the individuals in case group were DN cases; DM diagnosis and classification were according to World Health Organization (WHO) criteria, and DN was confirmed by the duration of DM, and the presence of urine albuminuria; evaluating the association between AGTR1 gene A1166C polymorphism and DN susceptibility; offering sufficient data on genotype distribution both in case and control groups; with reasonable grouping method, the case and control groups came from the same ethnical population, and were matched in gender and age; and focusing on human beings. Those publications were excluded from our study if they conformed to any one of the following conditions: case-only studies; with duplicated data; based on families or siblings; and letters, editorials, case reports, review articles, and conference abstracts.

Data extraction
Principal information of each eligible article was extracted independently by 2 reviewers, and contained first author's name, publication year, original country, ethnicity, genotyping method, numbers of cases and controls, genotype frequencies in case and control groups as well as P-value for Hardy-Weinberg equilibrium (HWE) in control group. If more than 1 study/ cohort were incorporated into 1 article, their data were extracted as separated ones. As for the disagreements on abstracted data, they were settled through discussion between the 2 reviewers; if no consensus was reached via discussion, a third reviewer would be consulted.

Quality assessment for eligible studies
Newcastle-Ottawa quality assessment scale (NOS score) was used to estimate the quality of the eligible studies. Eligible studies were classified into low, moderate, and high quality based on the NOS score 0 to 3, 4 to 6, and 7 to 9 scores.

Statistical analysis
All statistical analyses were completed with STATA 12.0 software (Stata Corporation, College Station, TX). The strength of the association between AGTR1 gene A1166C polymorphism and DN susceptibility was assessed by pooled odds ratios (ORs) with their corresponding 95% confidence intervals (95% CIs). Between-study heterogeneity was inspected with Chisquare-based Q test and I 2 test. When P < .05 and I 2 > 50%, significant heterogeneity presented, and random-effects model was used to calculate the pooled OR, otherwise, the fixed-effects model was used. Subgroup analysis based on ethnicity and DM type was performed to find the source of heterogeneity. Sensitivity analysis was conducted through sequential omitting each included study to test the stability of the final results. Begg funnel plot and Egger regression test were adopted to investigate the publication bias among the included studies visually and statistically, respectively.

Quantitative data synthesis
The main results of the meta-analysis are displayed in  (Fig. 2), and C versus A (OR = 1.208, 95% CI = 1.044-1.397) (Fig. 3) genetic models. Additionally, a similar influence of the polymorphism was also shown in Asian [under CC vs AA, CC + AC vs AA (Fig. 2)   AA contrasts] subgroups after stratification analysis by ethnicity and DM type. All these results illustrated that AGTR1 A1166C polymorphism was closely related to increased risk of DN in T2DM patients, especially in Asians.

Heterogeneity test
Q test revealed significant heterogeneity (P < .05) under all the 5 genetic comparisons, so the random-effects model was selected for calculating ORs. Then subgroup analysis based on ethnicity and DM type was performed to find the source of heterogeneity. There was no significant heterogeneity in the subgroups, indicating that the ethnicity and DM type might be the potential source of heterogeneity.

Sensitivity analysis
Sensitivity analysis was completed via deleting each selected study in turn to observe alteration in pooled ORs. During the whole process, no qualitative change occurred in the final results (data not shown), revealing the statistical robustness of our findings.

Publication bias examination
The shape of the funnel plots seemed symmetrical (CC vs AA + AC, Fig. 4), implying publication bias was negligible. Furthermore, these results were all confirmed by statistical data from Egger test (CC vs AA + AC, P = .142).

Discussion
DN is a common chronic complication, posting a great threat to healthy among DM cases. In recent years, the incidence rate of DN exhibits an upward trend in developed countries. In order to improve the management of DN, more and more researches are devoted to explore the pathogenesis of DN. Relevant data have confirmed that RAS possesses an important effect on the initiation and progression of DN via diverse mechanisms. Widely distributed among diverse tissues, AGTR1 mediates numerous vital biological effects, such as contracting vascular smooth muscle, stimulating the proliferation and thickening of vascular smooth muscle, and accelerating the release of aldosterone which are all related to blood pressure maintenance www.md-journal.com and target organ damages. As a G protein coupled receptor, AGTR1 has 7 highly conservative transmembrane functional domains. After binding with angiotensin, AGTR1 can activate G protein and 2 intracellular signal transduction pathways through inositol triphosphate and acetoglyceride. In one pathway, calcium releases activates protein kinases to promote the synthesis of proteins; while in the other pathway, cascade amplification of protein kinases activates MAPK which can prompt the expression of many protooncogenes after entering nucleus and thus further accelerating the division and proliferation of cells. In AGTR1 gene, there are 5 polymorphisms identified, namely A1166C, T573C, A1062G, G1517T, and A1878G, of which the former three are more common. The polymorphism A1166C is located at the 3 0 untranslated region of the gene, which has no effect on open reading frame and the encoding process of AGTR1 protein theoretically, but if it has linkage disequilibrium with adjacent chromosome positions with functional abnormalities, a series of impacts may occur in the stability of mRNA expression of the AGTR1 gene, AGTR1 number and distribution density as well as the affinity of AGTR1 to angiotensin II, thus strengthening the reactivity of angiotensin II and inducing the occurrence and development of DN. [18,19] Accumulating studies have discussed the relationship of the polymorphism A1166C with the susceptibility to DN, but no consistent opinion has been reached yet. For example, Doria et al [21] and Gallego et al [22] insisted there was no significant association between AGTR1 A1166C polymorphism and DN risk in Caucasian populations. However, the study by Shah et al [29] revealed a significantly much higher frequency of the C allele of the AGTR1 A1166C polymorphism in Indian DN patients, demonstrating the close relationship between the polymorphism and the disease in their studied population. In addition, Yin et al [20] in their research on Chinese also found that the C allele of the polymorphism is significantly more frequent in DN group than healthy control group and than DM without nephropathy group. All the discrepancies between the above findings might be partially attributed to different genetic backgrounds of participants in those studies, diverse selection criteria for study samples and uneven sample sizes. A system analysis was in urgent need to address the issues.
In order to obtain a reliable result about the genetic association of AGTR1 A1166C polymorphism with DN risk, the present metaanalysis was performed according to the guidance of PRISMA. After statistical analysis, the outcomes showed that AGTR1 A1166C polymorphism had significantly increasing-effect on DN susceptibility in total analysis, and a similar tendency was also revealed in Asian and T2DM groups after subgroup analyses of ethnicity and DM type. However, the significant association was not observed in Caucasian population and T1DM group. The results were partly consistent with the results of a similar meta-analysis carried out by Ding et al. [36] Their pooled analysis results demonstrated that AGTR1 A1166C polymorphism was obviously associated with DN in T2DM patients. Moreover, the significant association was not changed after stratification analysis by ethnicity. There were several reasons resulting in the divergences. Firstly, our meta-analysis included 17 eligible studies including 19 independent studies, while there were 10 studies included in the analysis of Ding et al. There were multiple recently published articles included in our analysis. Secondly, in our meta-analysis, patients in control group were all DM cases, however, 2 of the included studies in the pooled analysis of Ding et al set the healthy individuals as control. [37,38] The purpose of the pooled analysis was to investigate the genetic association of AGTR1 A1166C polymorphism with DN risk. The healthy individuals as control might cause bias to the final results. Thus, our meta-analysis included more high-quality articles that might provide more reliable and representative conclusions on these issues. In the present study, we investigated the effects of AGTR1 A1166C polymorphism on susceptibility of DN. The results obtained in our study might be helpful in identifying the population with high risk of DN among DM patients, especially among Asian and T2DM populations. However, these findings still need to be applied prudently due to several inevitable limitations in our study. The number of included studies was relatively small, which might result from source and language limitations in literature searching strategy. The relative small sample size might reduce the comprehensiveness of the final results. Moreover, possible combination and interaction of our studied polymorphism with other relevant factors were not embraced into the present study. Additionally, the complete medical records for the case and control populations were not available in all the included studies. The potential differences in clinical parameters might also cause bias to the final results. In view of the above mentioned restrictions in the present meta-analysis, these results need to be further verified in studies with larger sample sizes and more consideration of potentially collective effects.
In summary, our study displayed a risk-increasing influence of AGTR1 A1166C polymorphism on DN, especially among Asian and T2DM populations.


text:
Effect of Dietary Advanced Glycation end Products on Cell Cycle of Peripheral Blood Lymphocytes

Advanced glycation end products (AGEs) have been linked with autoimmune diseases. AGEs present in diet are known to be absorbed intestinally. We aimed to study effects of dietary AGEs present in repeatedly heated cooking oil on cell cycle profile of lymphocytes derived from healthy volunteers. AGEs were extracted by using Aqueous-TCAChloroform method and incubated with lymphocytes for 24 hours. All extracts showed ASF (15.42-75.18) between 355-440 nm wavelength. On incubation with AGEs, an increase in S-phase fraction of treated lymphocytes (7.67 ± 3.25) was observed in 10/10 samples analysed as compared to untreated lymphocytes (1.5 ± 1.6). The results suggest that dietary AGEs may have a role in triggering the proliferation of lymphocytes. Larger studies are required to explore role of dietary AGEs on properties of lymphocytes which may contribute to etiopathogenesis of autoimmune diseases. Such studies may provide us with a modifiable risk factor which may be altered to prevent autoimmune diseases. Citation: Chhabra A, Ram AK, Bhatia A, Goel S (2018) Effect of Dietary Advanced Glycation end Products on Cell Cycle of Peripheral Blood Lymphocytes. J Blood Lymph 8: 211. doi:10.4172/2165-7831.1000211

Introduction
Advanced glycation end products (AGEs) are the compounds formed due to non-enzymatic glycation of proteins, lipids or nucleic acids. Though, AGEs may form endogenously in the body as a result of normal metabolism, diet remains an important source of exogenous AGEs. Amongst the dietary AGEs, carboxymethyl lysine has been most commonly studied. Processing/cooking of foodstuffs under conditions of dry heat contributes to high AGE content of diets. Approximately 10-30% of AGEs present in the diet have been demonstrated to be absorbed intestinally. One third of absorbed AGEs are excreted via kidney and in stools. The accumulation of AGEs in body has been linked to several inflammatory and autoimmune disorders like rheumatoid arthritis, diabetes mellitus, myasthenia gravis, autoimmune uveitis etc. [1][2][3]. Although, indirect evidences suggest the role of dietary AGEs in etiology of inflammatory/autoimmune disorders a direct proof of principle for the same is lacking [4]. More recently, it has been established that consumption of AGE-rich diet leads to a significant rise in circulating AGE levels [5]. A reduction in dietary AGEs has also been shown to be associated with reduced oxidative stress and inflammation in patients with diabetes and renal disorders [6,7]. Most of the work on association of dietary AGEs and diseases has been done in diabetes and renal failure. In this work we have tried to study the effects of dietary AGEs on cell cycle profile of lymphocytes derived from healthy volunteers.

Materials and Methods
The AGEs were extracted from 10 samples of cooking oils subjected to repeated heating for frying different snacks by using Aqueous-TCA-Chloroform method. Briefly, 3 mL of oil sample was mixed and incubated with distilled water resulting in diffusion of water soluble AGEs in the aqueous layer. Subsequently, the proteins were precipitated out with Trichloroacetic acid (TCA, 13.5 M, 1 mL at 4°C for 10 minutes). The sample was then centrifuged (Thermo electron corporation; Sorvall®RC6PLUS) at 14,000 rpm for 10 minutes at room temperature. The supernatant was carefully separated from the protein pellet and equal volume of chloroform was added (to remove the lipids) followed by centrifugation at 14,000 rpm for 10 minutes at room temperature. The presence of AGEs in the extracts was confirmed by polyacrylamide Gel Electrophoresis and Per Iodic Acid Schiff staining which specifically stains glycoproteins and by determining AGE specific fluorescence (ASF) using fluorescence spectrophotometer (Elegant;Varian, Cary Eclipse, Excitation wavelength 355 nm and Emission wavelength 440 nm) [8]. The lymphocytes were isolated from 10 mL of blood samples by Ficolle Hypaque method. They were then incubated with AGEs (20 μL) extracted from cooking oil under culture conditions (37°C, 5% CO2) for 24 hours. At the end of the incubation period, effect on cell cycle was assessed by propidium iodide staining and flow cytometry (FACS Calibur San Jose USA) using CELL QUEST software. After evaluating at least 20,000 events, the percentage of cells in synthetic phase of cell cycle were taken as S-phase fraction.

Results
All extracts showed ASF at the wavelengths mentioned above. The ASF values ranged from 15.42-75.18 (Table 1). On incubation with AGEs, an increase in S-phase fraction (SPF) of the treated lymphocytes (Mean=7.67 ± 3.25) was observed in all 10 samples analysed as compared to untreated control lymphocyte population (1.5 ± 1.6) ( Figure 1). However, no relation was observed between ASF values and SPF.

Discussion
The effects of AGEs in autoimmune diseases have been attributed to activation of proinflammatory signalling pathways via the receptor for advanced glycation end products (RAGE). RAGE belongs to the immunoglobulin superfamily and is found on variety of cell types including macrophages, T-lymphocytes, endothelial, mesangial and smooth muscle cells. Besides AGEs, S-100 proteins, HMGB1, β-Amyloid and β-integrin also serve as ligands for RAGE. The AGE-RAGE axis most importantly activates NF-κB pathway leading to generation of proinflammatory cytokines like tumor necrosis factor-α, interleukin-6 and C-reactive protein [3]. Therefore, AGE-RAGE axis coupled with activation of downstream signalling pathways may play an important role in propagation of inflammatory responses. Infact, past studies conducted in animal models have demonstrated that blockade of RAGE helps to reduce inflammation in the joints involved by arthritis. RAGE has also been shown to lower the threshold for activation of TLR9 which is involved in recognition of self DNA immune complexes [4].
Our results suggest that dietary AGEs have a role in triggering the proliferation of lymphocytes. In a previous study on exogenous AGEs it was shown that binding of certain type of AGEs synthesized in laboratory [Bovine serum albumin (BSA)+D-glyceraldehyde and BSA+D-glycolaldehyde] to RAGE on monocytes can induce the expression of adhesion molecules, production of cytokines like IFN-γ, TNF-α and lead to proliferation of T-cells [9]. As FOXp3/RORγT balance is believed to be altered in autoimmunity, we believe that it will be interesting to explore the effect of dietary AGEs on it [10]. Since its discovery dysregulation of IL-33/ST2 axis has been implicated in pathogenesis of different inflammatory and autoimmune diseases. The above molecules are believed to be involved in differentiation of T helper (Th) cells to Th2, Th1 or Th17 depending upon the context and cell type. Investigating the influence of dietary AGEs on the above therefore appears to be a fascinating aspect in futuristic studies [11,12].
As dietary AGEs can be highly heterogeneous it remains to be determined which type of AGEs in diet are responsible for modulation of our immune response and whether such modulations can lead to autoimmune diseases. Further, which genetic, epigenetic and molecular events are exactly affected and what implications it may have in autoimmunity remains to be sorted out. Also, as AGEs constitute an  integral part of modern diets which include many heat processed food items, the threshold level beyond which they can be hazardous needs to be determined. Many studies have also shown the protective effect of substances like lime juice, vinegar etc. in diet. Marinating the foodstuffs with the later has been found to decrease the AGE content of foods.

Conclusion
This preliminary report suggests that dietary AGEs may have important effect on lymphocyte cell cycle and hence may serve as a basis for conducting larger scale studies to explore the role of dietary AGEs in etiology and promotion of autoimmune/inflammatory diseases. Such studies may be important as this may provide us with a modifiable factor which may be altered to curb the growing menace of autoimmune diseases across the globe.