Modeling of Mesoscale Variability in Biofilm Shear Behavior

Formation of bacterial colonies as biofilm on the surface/interface of various objects has the potential to impact not only human health and disease but also energy and environmental considerations. Biofilms can be regarded as soft materials, and comprehension of their shear response to external forces is a key element to the fundamental understanding. A mesoscale model has been presented in this article based on digitization of a biofilm microstructure. Its response under externally applied shear load is analyzed. Strain stiffening type behavior is readily observed under high strain loads due to the unfolding of chains within soft polymeric substrate. Sustained shear loading of the biofilm network results in strain localization along the diagonal direction. Rupture of the soft polymeric matrix can potentially reduce the intercellular interaction between the bacterial cells. Evolution of stiffness within the biofilm network under shear reveals two regimes: a) initial increase in stiffness due to strain stiffening of polymer matrix, and b) eventual reduction in stiffness because of tear in polymeric substrate.


Introduction
Bacteria can exist either as planktonic solutions or as surface associated colonies known as biofilms [1,2]. The latter is now recognized to be dominant mode of bacterial life. Biofilms generally include several differentiated populations of bacterial cells that are embedded in a matrix of self-produced extracellular polymeric substances (EPS). EPS acts as natural 'glue' providing mechanical integrity to the biofilm structure [3]. Due to this social form of growth biofilm bacteria offers resilience to external stresses and thus provides itself with an existential advantage over planktonic form of living [4][5][6][7]. Biofilms have been implicated for their role in recurrent infections and biofouling [8,9]. However, they also play a crucial role in several environmental processes including chemical cycles and waste-water treatment [10,11]. Our need to control biofilm growth, as well as possess the ability to remove undesirable biofilms, necessitates a proper understanding of the material characteristics of this biological soft matter [12].
Even though there exist other polymeric materials within the EPS matrix that shows shearthinning behavior for small time scales, only the shear-stiffening protein chains have been considered in the present study as a first approximation [40]. Since the bacteria located within the biofilm network acts as a rigid body, with increasing bacteria loading, effective stiffness of the biofilm network increases. The assumption of considering bacteria as a rigid object is based on prior experimental works by other researchers [41,42]. In the present context, "bacteria loading" indicates the volume fraction of the biofilm domain occupied by living bacteria. In fact in the context of soft matter physics, rigid bacteria are often considered analogous to non-deformable colloidal particles [6]. While this assumption might not apply to all possible scenarios, we consider the domain space where this assumption is valid. Different portions of the same biofilm can have various bacteria loading. Fig 1(b), 1(c) and 1(d) shows three different cases of 26%, 42% and 55% bacteria loading, respectively. These experimentally observed images of biofilm belong to Escherichia Coli 042 bacteria grown in the presence of an electric field, in a miniature culture chamber [43].
Under shear deformation, the stress-strain curve obtained for different bacterial loading has been shown in Fig 1(e). With increasing proportion of rigid bacteria, stiffness of the system increases. Strain stiffening is observed under all the three bacterial loading conditions. The spring stiffness parameter has been assigned in such a way that the stress-strain curve for the 26% bacteria loading corresponds to the experimental result reported by Stoodley et al. (2002) [5]. The exact parameters used for this comparison are reported in Table 1. Under externally applied shear load, the biofilms behave as elastic solid for a short period of time, and shows viscoelastic response at longer time intervals [5]. In the present context, only elastic response of the biofilm is being simulated. The stress-strain data-points adopted from Stoodley et al. were obtained by changing the externally applied load at a time interval of 5 seconds. The mechanical response at such short time intervals (with respect to the viscoelastic relaxation time scale of a biofilm) can be assumed to be elastic [13]. These biofilms were generated from Pseudomonas aeruginosa strains and grown under laminar flow conditions for 14 days with sufficient nutrient feed. No specific information was provided whether the biofilm being tested were at the beginning or end of life. According to the analysis reported in Schaudinn et al. [44], biofilms containing Staphylococcus Epidermis cultivated in static liquid cultures, that is 14 days old should be considered to have reached its end of life. The EPS matrix usually consists of polysaccharides, proteins, nucleic acids, phospholipids and other humic substances [45]. In the present context, it is assumed that the EPS matrix consists of polymeric protein chains that can uncoil under externally applied load.
It can be observed from Fig 1(e) that for higher bacteria loading, the strain stiffening starts earlier. An important point to note is that in our model, bacteria have been approximated as rigid particles [6]. Hence, under the application of an external force, majority of the deformation occurs within the EPS matrix because of its smaller stiffness. For the sections of biofilm with high bacteria loading also corresponds to low EPS content. All the strain stiffening comes only from the EPS matrix. To maintain equilibrium, all the portions of the biofilm should experience equal force. Due to the rigidity of bacteria, small deformation leads to large force. To produce similar amount of force in the EPS matrix, significant deformation has to occur. Hence, the springs adjacent to the bacteria enters the strain stiffening regime early enough. For biofilms containing large amount of EPS matrix, the strain gets equally distributed among all the protein chains, and strain stiffening initiates late in each of them. Hence, the biofilm microstructure that contains smaller proportion of EPS, experiences more strain, and the stiffening initiates earlier. On the other hand, polymeric substrates in biofilm microstructures containing a higher proportion of EPS are subjected to less strain, leading to delayed initiation of strain stiffening phenomena. (e) Shear stress-strain response of the biofilm for the three types of bacteria loading scenario. Since the bacteria are considered as rigid objects, larger volume fraction of bacteria leads to higher stiffness for the biofilm. Stiffness for each lattice spring element in the EPS matrix has been set in such a way that the biofilm with 26% bacteria loading follow the experimental shear stress-strain data obtained from Stoodley et al (2002). (f) Based on the stiffness of the EPS matrix and the bacteria, Hashin-Shtrikman (H-S) lower and upper bounds of shear stiffness for the biofilm has been evaluated. Shear stiffness of the biofilm at 5% and 35% shear strain falls within the Hashin-Shtrikman bounds. The parameters used in these simulations are given in Table 1. To check whether the stiffness values estimated using DBM is correct or not, a bound based analysis has been conducted in Fig 1(f). The Hashin-Shtrikman (H-S) upper and lower bounds for shear modulus were plotted along with stiffness of the biofilm for all the three different bacteria loading at 5% and 35% shear strain [46]. The exact expressions of Hashin-Shtrikman bounds have been estimated by using a variational approach and applicable to any random geometry of the bacteria inclusions. The Hashin-Shtrikman bound has been used to verify the biofilm elastic modulus predicted by DBM. The lower and upper bounds of the shear modulus are denoted by G l and G u , respectively. Their expressions are given as follows: and Here, G represents shear modulus, K represents bulk modulus and c indicates the volume fraction. The subscript "1" and "2" signifies properties of the EPS matrix and the bacteria, respectively. The stiffness at 5% strain corresponds to the lower bound, whereas; at 35% strain the shear stiffness lie perfectly inside the H-S bounds. This implies that the H-S bounds for shear stiffness are applicable in biofilm networks as well for the purpose of checking the validity of experimental results. These H-S bounds are applicable to both 2D and 3D two-phase media. Since the 2D biofilm geometries simulated here are plane strain approximations of a 3D geometry with zero stress along the out-of-plane direction, the H-S bounds are applicable. In the present 2D analysis, it is assumed that similar volume fraction and morphology of bacteria is observed along the out-of-plane direction of the biofilm. Under these conditions, the 2D analysis reported here is capable of capturing the 3D response of the realistic biofilms. However, if the volume fraction or morphology of bacteria changes along the out-of-plane direction, approximation of the entire 3D biofilm by a 2D representation is not possible.
It is conjectured that strain stiffening occurs within biofilm networks because of the unfolding of protein chains. The question is whether the entire EPS matrix experiences uniformly distributed shear strain, or localization of shear strain occurs within the DBM network. The biofilm morphology taken into consideration consists of 42% bacterial loading with uniform distribution (also described in the methodology section). Fig 2(a) shows (in cyan) the location of those springs that has experienced unfolding. It has been observed that very small amount of unfolding occurs on the left and right side of the biofilm network. Majority of the strain got localized along the diagonal direction, which led to unfolding of springs along the diagonal. This type of shear strain localization along the diagonal has also been observed in metals [47]. From this analysis, it is evident that unfolding of the protein chains is localized along the

Name of parameter Unit Value in EPS matrix Value in rigid bacteria
Ref.
Young's modulus N/m 2 1.0 100.0 [5] Poisson's ratio -0.33 0.33 -Shear modulus N/m 2 0.33 33.33 [5] Average energy threshold for strain-hardening J/m 2 10.0 - [57,58] Energy threshold for sub-spring uncoiling (applicable to each spring) J 1.25e-5±0.25e-5 -- diagonal of the EPS matrix. Rupture behavior has not been simulated in this case. The parameters adopted to run the simulations provided in Fig 2(a), 2(b) and 2(c) are provided in Table 1. Force contours observed inside the biofilm network have been plotted in Fig 2(b) and 2(c) under externally applied strains of 5% and 35%, respectively. For 5% strain, small amount of force has been observed within the network, which has a maximum value of 0.3N. As the strain is increased to 35%, large amount of force is generated, which shows a maximum value of 4.5N. It is clear that seven fold increment in strain (5% to 35%) results in fifteen-fold increase in the magnitude of force (0.3N to 4.5N). This nonlinearity is attributed to the strain stiffening behavior observed due to the unfolding of EPS protein chains. Since unfolding of sub-springs lead to stiffening of the springs, an unfolded spring produces larger force. Thus, the unfolded springs along the diagonal direction observed in Fig 2(a) naturally leads to localization of force along the diagonal (Fig 2(b) and 2(c)). Force contour observed at 35% shear strain. Magnitude of force is around 4.5N, which is an order of magnitude greater than the forces observed under 5% strain. Larger amount of forces accumulate along diagonal directions, which eventually leads to enhanced strain hardening along the diagonal direction. Similar behavior has been reported in part (b) of this image. The parameters used in these simulations are given in Table 1.
None of the simulations discussed till now had rupture of the springs incorporated within the DBM. The biofilm microorganisms (such as, bacteria) are embedded within the EPS matrix which is composed mostly of polysaccharides, proteins and lipids [48]. As rupture of springs is added into the simulations, the EPS matrix may experience detachment from the bacteria. The EPS matrix interacts with the environment by attaching biofilms to a surface. Sorption and transport of dissolved and particulate substances from the environment provides nutrients for biofilm organisms [44,45,48]. Rupture of the EPS matrix can affect the chemical communication within the bacteria, which has the potential to impact the lifespan of a biofilm. In a nutrient poor condition, the bacteria may escape to the planktonic phase. By incorporating rupture of the protein chains, estimation will be made about the conditions when tearing can appear inside the EPS matrix. To simulate various bacterial morphologies, we generated various bacterial flocculation based on experimental images [43], which were later digitized. We considered various scenarios, such as, bacteria can be uniformly distributed (as shown in  Table 2 have been adopted. Density of ruptured protein chains within the EPS matrix is relatively less for the uniform distribution (see Fig 3(g)). As a result some portion of every bacteria are connected with the EPS matrix. On the other hand, the biofilm morphology with centrally clustered bacteria experiences larger amounts of rupture inside the EPS matrix adjacent to the bacteria (see Fig 3(h)). Extreme amount of tear within EPS has the potential to completely detach some of the bacteria from the biofilm network, which can impact their chemical communication. The amount of rupture and tear in the EPS matrix within uniformly distributed and clustered biofilm seems minimal. However a thorough analysis clearly indicates that several bacteria within the clustered sample (shown in Fig 3(h)) is completely detached from the EPS matrix and suffers due to lack of nutrients. All the bacteria are connected to the EPS matrix for the uniform distribution (shown in Fig 3(g)). Thus, high clustering of bacteria can be beneficial for removal of biofilm through the application of mechanical force. The physical reason behind observation of enhanced rupture in high clustering, stems from the fact that regions with high bacterial loading also experience lower concentration of EPS. Majority of the externally applied strain is accommodated by the soft EPS matrix. Strain stiffening is observed earlier in regions with less amount of polymeric substrate, which subsequently leads to rupture. This phenomenon is similar to that observed in Fig 1(e).
Due to the presence of larger volume fraction of EPS matrix adjacent to each bacterium, less strain stiffening and rupture occurs under uniform distribution. The corner clustering case has been investigated to understand the propagation of tear through the EPS matrix in between two clusters. It is evident from Fig 3(i) that the cracks through the EPS matrix propagates from one cluster to the nearest neighboring cluster. Probability of propagation of tear along the diagonal direction is relatively less. Rupture propagation is a mechanism to release strain energy. Connecting the ruptures in the nearest neighboring clusters provides a shorter pathway to release strain energy quickly. That is why no rupture patterns have been observed which connects diagonally opposite clusters.
Stress vs. strain response of the biofilm network under shear strain induced loading has been shown in Fig 4(a) using the parameters reported in Table 2. Mechanical rupture has also been incorporated within the simulation. Due to unfolding of springs shear stiffness increases initially. Once the rupture of EPS matrix is initiated, stiffness of the network decreases significantly. Stress-strain curve for both the central clustering and uniformly distributed bacteria loading has been shown. The result obtained from computational analysis corresponds very closely with that observed in experiments conducted by Korstgens et al. (see Fig 4(a)) [19]. Because of localized strain stiffening, early initiation of rupture has been observed in the central cluster. For this particular case, cracks initiate at a stress level of 1100N/m 2 . For uniformly distributed bacteria loading, rupture initiation starts at a stress level of 1300N/m 2 . Delay in microcrack nucleation happens because of the fact that uniformly distributed bacteria contains large amount of EPS matrix adjacent to it, which acts as a buffer against rupture initiation. (g) Distribution of rupture within the uniformly distributed bacteria. When all the bacteria are well dispersed, damage in the biofilm is not localized enough for an entire bacterium to be completely detached from the EPS matrix. (h) Rupture profile as observed for centrally clustered bacteria. When the bacteria are clustered within the biofilm, force concentration is very high around the bacterium, which leads to mechanical detachment of the EPS matrix from certain bacteria. (i) Evolution of damage for corner clustering scenario. Purpose of the corner clustering is to investigate how the tear would propagate through the EPS between two clusters of bacteria. It can be concluded that cracks through EPS matrix propagates to the nearest cluster. There is very little possibility for a tear through EPS to span between two clusters which are located far away from each other. The parameters are given in Table 2. Evolution of shear stiffness has been plotted in Fig 4(b) with respect to increasing shear strain. It has been observed that strain stiffening dominates till a particular strain. Because of the presence of tear in the protein chains of EPS matrix, rupture behavior starts to dominate after that critical strain. Two regions where strain stiffening dominates, and where rupture dominates have been clearly marked in Fig 4(b). For the case of central clustering, domination of rupture over strain stiffening starts earlier than the case of uniformly distributed bacteria. This early initiation of fracture in the EPS matrix is attributed to the fact that in centrally clustered biofilm, very little amount of EPS matrix exists adjacent to the bacteria. Under shearinduced deformation, the EPS matrix carries majority of the strain and experiences early initiation of strain stiffening and subsequent rupture in the regions with less amount of the polymeric substrate. As a result, the uniformly distributed bacteria experiences delayed initiation of strain stiffening and rupture due to the presence of excessive amount of EPS matrix around each bacterium.   Table 2. doi:10.1371/journal.pone.0165593.g004

Discussion
In our effort to create the DBM framework, we have made certain simplifying assumptions. For example, the biofilm has been assumed to be a two phase system, soft EPS matrix and rigid bacteria. However, in a realistic situation EPS matrix can contain inert biomass, water and/or air filled voids [5,19,[49][50][51]. As a result, practical situations are quite complex and a full model has to incorporate these different aspects. Void spaces filled with air would be of paramount interest in biofilms that are partially exposed to air (such as those grown on Agar plates). In the present context, images that motivated the DBM model were taken from biofilms that are totally immersed in a liquid environment and hence air-like voids have not been considered [43]. It is assumed that the biofilm had sufficient time to equilibrate with its aqueous environment and hence the EPS/water network structure could be considered to be a homogenous system. Extending the current computational methodology to such complex situations is indeed possible, but would require extensive validation and is left as an exercise for a future work. Morphology of the bacteria cells significantly impact the overall viscoelastic and dynamic behavior of biofilms [37]. Another simplifying assumption worth mentioning is the fact that effect of bacterial cell morphology has been neglected in the present study. However, the shape of the cells adopted here has been inspired by experimentally observed images of biofilms [43]. Since the protein chains within the EPS matrix mostly govern the strain-stiffening behavior, bacterial morphology has minimum impact there. Hence the effect of cell morphology can be neglected as a first approximation. Detailed analysis of the shape of the cells will be conducted in future.
While it is known that biofilms can exhibit complex viscoelastic response to external forces, only elastic response has been modeled here. This is equivalent to considering a situation where the response time (t) is much smaller than the viscoelastic relaxation time-scale (λ). Some studies indicate that biofilms can be characterized by a dominant viscoelastic time scale of 20 minutes [52] and hence considering a scenario where t << λ as elastic, is quite realistic. Moreover, viscoelastic responses are often linked to very large deformations [53], which is outside of the scope of the current investigation.

Conclusion
In this article, a novel digital biofilm model has been developed, which has the potential to take into consideration the effect of any bacteria-EPS microstructure observed in real biofilms. Different volume fractions of bacterial cells show that for higher bacteria loading, stiffer stressstrain response is obtained. Also, under high bacteria loading enhanced strain hardening occurs within EPS. It leads to earlier initiation of tears inside the protein chains. Similar to strain localization observed in metals, localized strain stiffening along the diagonal direction has been observed in biofilm networks under shear deformation. Enhanced rupture has been observed in EPS matrix that is adjacent to the bacteria in a closely packed cluster. As a result, detachment of the bacterial cells from the EPS matrix is possible in a highly clustered domain under mechanically induced shear type load. By plotting the evolution of stiffness in a biofilm network, two different regions have been recognized. In the first domain strain stiffening dominates and the effective stiffness of the biofilm network increases. In the second region, the effect of rupture of the polymer matrix plays a major role and eventually results in reduction of the effective stiffness. As a living system biofilm adapts with time and its biological activity is intricately tied to its mechanical and transport behavior. Most of these are currently neglected as we have considered only the short-time scale elastic strain-hardening and rupture response of biofilms to external shear. This implies that the time-scale of shear (τ) is much smaller than the relaxation time scale of the biofilm (λ), i.e. τ << λ. It has been reported that λ~20 minutes for many different biofilm species [52]. Hence for shear forces lasting a short time, it is an appropriate first approximation to neglect growth and other biological activities [13]. Elastic rupture of the EPS matrix in the shorter time scale allows for the viscoelastic movement of biofilms in the longer time range. Naturally, for simulating much longer time scales the assumption of linear elastic behavior would not remain valid.

Methodology
In this article, morphology aware computational model has been developed that can elucidate the shear behavior of a biofilm. Here the word "morphology aware" indicates the fact that the developed model can simulate any morphology of the bacteria cells as well as different clustering of bacteria within the biofilm (such as, uniformly distributed or clustered presence of bacteria cells). In order to investigate the effect of shear force on a biofilm (Fig 1(a)), we begin by using experimental visualization of the biofilm structure to create a digital representation of the biofilm (Fig 1(b)-1(d)). Once the digital representation of the biofilm is complete, then a single domain modeling approach is applied by superimposing the digitized biofilm microstructure on top of a network of lattice springs (see Fig 5(b)). Here the word "single modeling approach" indicates that only a single mesh has been generated. Different regions of the mesh has been assigned different properties based on the whether it belongs to a "bacteria" or "EPS matrix". This computational characterization of a realistic biofilm microstructure has been named as Digital Biofilm Model (DBM).
As shown in Fig 5(a), in the digital biofilm model a triangular network of lattice springs has been considered to analyze the mechanical deformation of the biofilm. In the quasistatic computational scheme adopted in the current research, entire mass of the system is assumed to be concentrated at uniformly distributed nodes [39]. Any two adjacent nodes are connected by a spring, which shows only axial resistance. The triangular spring network leads to a coordination number of six. Interaction with only the nearest neighbors has been taken into consideration. Because of the triangular configuration, Poisson's ratio of 1/3 has been observed for this network.
Superimposing the digitized biofilm on top of the triangular network, springs that lie within the bacteria and springs that reside inside the EPS matrix has been separated (see Fig 5(b)). In order to estimate stiffness of the biofilm network, the bacteria were assumed to be rigid with respect to the extracellular polymeric substrate (EPS) [6]. Furthermore, we restrict ourselves to studying biofilm response to shear forces for time scales much smaller than the viscoelastic relaxation time scales (i.e. t << λ) [5]. The elastic behavior dominates biofilm response in such small time scale. The EPS matrix consists of protein chains, which show strain-stiffening behavior due to unfolding of the polymer-chains [30]. Since the biofilm consists of soft polymeric gels, due to external effects the film can go through large deformations. Sometime shear strains as large as 50% can be observed because of flow induced force [7].
The rules applied to the spring are: (a) springs which are representative of bacteria have a very high spring constant [6], (b) springs which are in the EPS have a low spring constant and an energy prescribed limit is imposed. When the threshold strain energy is exceeded the springs become 'uncoiled' and are represented by a higher magnitude of spring constant. This latter rule is inspired by findings of Buehler [29,30]. Using molecular dynamic simulations, unfolding of the helical patterns were captured which leads to stiffening of the entire polymer chain. The untangling mechanism was initiated by the stutter defect present within a protein chain. Once all the protein chains go through the unfolding mechanism, the extended fiber show relatively high stiffness, but cannot sustain a significant rise in externally applied strain. (c) Once the strain energy stored in a spring exceeds its rupture threshold, the spring breaks. It is irreversibly removed from the network of lattice springs and the force carried by this spring gets distributed among its neighbors [54].
To capture the unfolding of springs, it is assumed that each spring in the network consists of eight folded sub-springs [32,55,56]. Due to application of shear strain, as each of the subsprings unfold stiffness of the full spring increases (see Fig 5(c)). After the eighth sub-spring has unfolded, application of more external load leads to rupture of the entire spring [38]. The remaining question is how to determine when the unfolding of sub-springs will occur. An energy prescribed criteria has been adopted in this study. Under externally applied loads, the local force (f) and displacement (u) has been calculated which equilibrate the spring network. Strain energy in each spring is calculated as, ψ = (1/2)fÁu. With increasing amount of externally applied load, the total strain energy stored in each spring increases. Eight uniformly distributed random energy thresholds have been generated for each spring (ψ t,i ;I = 1,2,3,. . ., 8) in such a way that ψ 1 < ψ 2 < ψ 3 <. . .< ψ 8 . As the strain energy in the full spring exceeds an unfolding threshold (ψ > ψ t,i ), one sub-spring unfolds and the stiffness of the full spring is increased (shown in Fig 5(c)). This untangling mechanism continues until all the eight sub-springs have unfolded. Stiffness of the spring reaches its maximum value at that point [30].
Another energy threshold has been assigned to all the springs that correspond to its ultimate strength (ψ t,u ) and it follows the constraint ψ t,u > ψ 8 . Once the strain energy in the completely unfolded spring exceeds its ultimate rupture threshold (ψ > ψ t,u ), the spring is assumed to be broken. Hence, it is irreversibly removed from the lattice network and the contribution of the spring is completely eliminated from the stiffness matrix. Rupture in the lattice spring network results in reduction in stiffness of the entire biofilm. Since the spring constant of the full spring increases with each unfolding of the sub-springs, the stiffness of the lattice network also changes. In order to correctly capture the increase in stiffness and model the spring deformation accordingly, an iterative technique has been adopted. The externally applied load is increased in increments of very small magnitude. All the sub-springs that unfold or full springs that break within a loading interval are modified accordingly. The energy in each springs are also estimated in increments, ψ n+1 = ψ n + (1/2)ΔfÁΔu. Here, ψ n+1 and ψ n signifies the strain energy at the current and the previous steps, respectively. Δf and Δu corresponds to the incremental force and displacement, respectively, observed in the current loading interval.
Estimation of the correct parameters of the EPS matrix and the rigid bacteria that can successfully predict the strain hardening and rupture behavior of the biofilms is a challenging task. However, the elastic behavior, in terms of stiffness, strain hardening and rupture, of an entire biofilm has been reported in some of the recently published articles [5,6,19]. For the simulations reported in this article, the elastic properties of the EPS matrix and rigid bacteria have been estimated such that the overall biofilm response closely follows the experimentally observed pattern. Exact values of the parameters adopted for the comparison with strain hardening and rupture behavior of biofilms are provided in Tables 1 and 2, respectively. The main purpose of this article is to demonstrate the basic physics behind the strain hardening and subsequent softening observed in biofilms. The range of parameters adopted to study these effects does not alter the physical explanation. Hence, usage of highly varying magnitude of parameters is well justified. Detailed discussion regarding the mathematics behind the computational methodology and a study of mesh size dependence has been incorporated within the Supplementary section (see S1 File).