A study of nitric oxide dynamics in a growing biofilm using a density dependent reaction-diffusion model

A study of Nitric Oxide dynamics in a growing biofilm using a density dependent reaction-diffusion model

Journal Pre-proof

A study of Nitric Oxide dynamics in a growing biofilm using a density dependent reaction-diffusion model Maryam Ghasemi, Benjamin Jenkins, Andrew C. Doxey, Sivabal Sivaloganathan PII: DOI: Reference:

S0022-5193(19)30422-9 https://doi.org/10.1016/j.jtbi.2019.110053 YJTBI 110053

To appear in:

Journal of Theoretical Biology

Received date: Revised date: Accepted date:

25 December 2018 5 July 2019 15 October 2019

Please cite this article as: Maryam Ghasemi, Benjamin Jenkins, Andrew C. Doxey, Sivabal Sivaloganathan, A study of Nitric Oxide dynamics in a growing biofilm using a density dependent reaction-diffusion model, Journal of Theoretical Biology (2019), doi: https://doi.org/10.1016/j.jtbi.2019.110053

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

1

2

3 4 5

6 7 8

9 10 11 12

13 14 15 16 17

18 19 20

Highlights • • Oscillatory behavior of N O• concentration under a microaerobic regime is a temporary phenomenon and does not generate gradients in N O• concentration levels within the biofilm. • In order to weaken the defense system against N O• (i.e. maximize the inhibitory effects of N O• against bacterial biofilms), large amounts of N O• donor with short half-lives should be used to treat the biofilm. • The initial size of a biofilm colony affects the activity of N O• reductants. By decreasing the initial size of biofilm colonies, less reducing agent is produced hence more N O• remains in the system which consequently removes biofilms effrciently. • Nutrient deprivation stress can change the N O• dynamics depending on the oxygen and N O• concentration and biofilm size. Thin biofilms respond to starvation stress under microaerobic conditions. In this case, defense system against N O• becomes weak and more biofilms can be eradicated by N O•. • For thick biofilms under microaerobic regimes, starvation prevents N O• from diffusing into the inner region of the biofilms and leads to the layering phenomenon.

1

21

22

23

24

A study of Nitric Oxide dynamics in a growing biofilm using a density dependent reaction-diffusion model Maryam Ghasemi1,a , Benjamin Jenkins2 Andrew C.Doxey2 , Sivabal Sivaloganathan 1

1

Dept. of Applied Mathematics,

Univ. Waterloo, Waterloo, ON, Canada, N2L 3G1 2

Dept. of Biology,

Univ. Waterloo, Waterloo, ON, Canada, N2L 3G1 a

25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43

corresponding author

Abstract One of a number of critical roles played by N O• as a chemical weapon (generated by the immune system) is to neutralize pathogens. However, the virulence of pathogens depends on the production activity of reductants to detoxify N O•. Broad reactivity of N O• makes it complicated to predict the fate of N O• inside bacteria and its effects on the treatment of any infection. Here, we present a mathematical model of biofilm response to N O•, as a stressor. The model is comprised of a PDE system of highly nonlinear reaction-diffusion equations that we study in computer simulations to determine the positive and negative effects of key parameters on bacterial defenses against N O•. From the reported results, we conjecture that the oscillatory behavior of N O• under a microaerobic regime is a temporal phenomenon and does not give rise to a spatial pattern. It is also shown computationally that decreasing the initial size of the biofilm colony negatively impacts the functionality of reducing agents that deactivate N O•. Whereas nutrient deprivation results in the development of biofilms with heterogeneous structure, its effect on the activity of N O• reductants depends on the oxygen availability, biofilm size, and the amount of N O•.

2

44

keywords: Nitric oxide, Biofilm, Nonlinear model, Simulation, Oscillation

45

Mathematics Subject Classification: 92D25, 35K65, 65M08, 35K59

46

1

Introduction

47

The emergence of resistance of bacteria to antibiotics has intensified the

48

search for alternatives to fight infections which do not inadvertently promote

49

bacterial growth and increase pressure on the host environment [40]. One

50

potential class is antivirulence therapy which targets the host-pathogen in-

51

teractions that are vital to pathogenesis [10]. Inhibiting the determinants

52

of virulence such as biofilm formation, quorum sensing and detoxifying the

53

immune system generated antibacterials limits the pressure to the host en-

54

vironment which results in less resistance compared to the application of

55

conventional antibiotics [4]. N O•, a potent antimicrobial agent, is a highly

56

diffusible gas generated by mammalian N O• synthesis (NOS) to neutralize

57

pathogens [2,47]. The bacterial defense system against N O• has been identi-

58

fied as a virulence factor [39]. In order to impair the bacterial defense system

59

and find appropriate approaches to inhibit N O• detoxification, identifica-

60

tion of the components that react with N O• and an understanding of the

61

underlying kinetic process involved are essential.

62

Although identifying inhibitors with more favourable properties than the

63

existing known inhibitors can be acheived by performing chemical screen,

64

the study of the bacterial N O• response mechanism (to recognize the other

65

components that can be targeted to enhance the susceptibility of bacteria) 3

66

is complicated. This is because of the multifaceted role of the N O• stress

67

response mechanism and the broad reactivity of N O• (including the reac-

68

tion with iron-sulfur ([F e − S]), DNA and thiols) [2, 21, 26, 35, 40]. Exposing

69

bacteria to N O•, stimulates different types of activities inside the bacte-

70

ria. Some biomolecules are damaged while others are repaired and various

71

other activities are inhibited. This gives a complex reaction network de-

72

termining the outcome of N O• and its effect on the bacterial culture. For

73

quantitative study of this network, various kinetic models have been pro-

74

posed at different levels of complexity to predict the outcome of N O• treat-

75

ments [19, 21, 25, 26, 37, 38, 41, 42]. However, in general, they have considered

76

only the temporal behavior of N O• inside a well mixed bacterial culture.

77

Besides the deficiency of kinetic models in describing N O• dynamics in cul-

78

tures with nutrient gradients and spatial distribution of N O•, the models

79

are large systems of Ordinary Differential Equations (ODEs) with many pa-

80

rameters, with tenuous links to the actual physical problem. This makes it

81

difficult and cumbersome to identify the key parameters that influence N O•

82

detoxification.

83

It is well known that bacterial cultures do not always grow as planktonic

84

cells. In fact, they can also emerge in clusters known as biofilms. Bacterial

85

biofilms are microbial communities attached to an immersed surface known as

86

the substratum. A main feature of biofilms is that they produce gel like layers

87

known as Extracellular Polymeric Substances (EPS) in which bacteria are

88

embedded and which offers them protection against chemical and mechanical 4

89

washout and against antimicrobial agents [18, 22, 49].

90

To understand the process of biofilm growth and to predict its behavior

91

under different conditions, several mathematical models have been proposed

92

in the literature. An early prototype biofilm model proposed in [6] (which

93

has been derived both from the viewpoint of biofilms as spatially structured

94

populations [20] and as a fluid [12, 32]), is perhaps the first model which

95

reflects the ecological-mechanical duality of biofilms. In its basic form, as a

96

biofilm growth model, it is a nonlinear density-dependent reaction-diffusion

97

equation for biomass that is coupled with a semi-linear reaction-diffusion

98

equation for growth limiting nutrients. The biomass equation shows two in-

99

teracting nonlinear diffusion effects: (i) porous medium degeneracy as the

100

dependent variable vanishes, and (ii) a super-diffusion singularity as the de-

101

pendent variable in the diffusion coefficient approaches the known maximum

102

density. The interplay of these two nonlinear effects ensures that (1) the

103

solution of the biomass equation is bounded by the maximum cell density re-

104

gardless of growth activity [9], (2) that spatial expansion of the biofilm does

105

not take place if there is enough space for new biomass to be accumulated,

106

and (3) that a distinct interface between the biofilm colonies and the sur-

107

rounding aqueous phase is established. Since biofilm colonies are subject to

108

growth, these interfaces are not stationary but change over time. Eventually

109

neighbouring colonies may merge into a bigger colony, in which case their

110

interfaces can merge and dissolve as well.

111

Whereas mathematical models have been proposed to study the N O• 5

112

stress response mechanism in planktonic cultures, there has been no cor-

113

responding attempts in the literature to investigate the dynamics of N O•

114

within the biofilms, primarily due to the complexity of their structure. Our

115

main objective is to focus on this latter problem, and for this purpose, we

116

refine the prototypical single-species biofilm model to account for the ma-

117

jor components in the N O• detoxification system (under aerobic, anaerobic

118

and microaerobic conditions). In our model, the biofilm consists of species

119

that are major participants in the defense system against N O•. Dissolved

120

substrates are glucose as the only carbon source, oxygen and an N O• donor

121

such as dipropylenetriamine (DPTA) that generates N O• exogenously. The

122

resulting density dependent reaction-diffusion model is too complex to per-

123

mit rigorous mathematical analysis of qualitative model behavior and so we

124

probe the model through computer simulations to try and predict N O• dy-

125

namics in the growing biofilm. This helps identify the key parameters that

126

play a role in inhibiting determinants of virulence and thus contribute to

127

increasing the susceptibility of the biofilm to N O• treatment.

128

2

Mathematical Model

129

2.1

Assumptions to be included in the model

130

We develop a quasilinear spatio-temporal model that describes the dynamics

131

of N O• in a growing biofilm. For this purpose, the following assumptions

132

are made: 6

133

1. The biofilm is defined as a region with positive biomass density which

134

is surrounded by an aqueous phase in which the biomass density is

135

zero [6,7]. The biofilm consumes the nutrient and expands. We assume

136

that nutrient is dissolved and transported in the liquid and biofilm

137

phase by Fickian diffusion.

138

2. We assume that the biofilm, more specifically Escherichia coli (E.coli)

139

biofilm, consists of five species which are N O• reductants: (1) flavo-

140

hemoglobin (Hmp) which converts N O• to nitrate under aerobic con-

141

ditions and protects bacteria against N O• toxicity; (2) reduced fla-

142

vorubredoxin (N orV ) that converts N O• to nitrous oxide (N2 O) and

143

water under anaerobic conditions; (3) cytochrome bd − I, and bd − II

144

(referred to collectively as Cyd) as respiratory oxidase. Respiratory ox-

145

idases are the last enzymes in the aerobic respiratory chain and found in

146

many bacterial pathogens. Expression of Cyd is a defence mechanism

147

in response to N O• under microaerobic conditions with the oxygen con-

148

centration levels between 5µM to 35µM [45]. (4) complex Cyd − N O•

149

which is produced after reversible binding of N O• to Cyd; (5) biofilm.

150

Adequate nutrient supply results in biofilm growth which consequently

151

increases biomass volume fractions (i.e. N O• reducing agents) depend-

152

ing on the oxygen concentration and upon the existence of N O• (which

153

acts as a stressor and stimulates the growth of biomass volume frac-

154

tions). 7

155

3. It is assumed that nitric oxide is generated exogenously by treating the

156

biofilm with DPTA NONOate (dipropylenetriamine NONOate) (Z) −

157

1−[N − (3 − aminopropyl) − N − (3 − ammoniopropyl)amino] diazen−

158

1−ium−1, 2−diolate as N O• donor. DPTA dissociates with a half-life

159

of ∼ 2.5 h (at 37◦ C and pH 7.4) to release 2 mol of N O• per mol of

160

parent compound. Moreover, it is assumed that N O• kills the bacterial

161

biofilm but is not consumed by biofilm.

162

4. Dissolved substrates N O•, oxygen and DPTA are transported by dif-

163

fusion in the surrounding liquid and in the biofilm through Fickian

164

diffusion. For the sake of simplicity, we assume that diffusion coeffi-

165

cient of these substrates is not density dependent and that they diffuse

166

homogeneously.

167

5. It is also assumed that Hmp and N orV decay naturally.

168

2.2

Governing Equations

169

Based on the above assumptions, the mathematical model is formulated as

170

a system of nonlinear PDEs over a domain Ω ⊂ R2 . The dependent vari-

171

ables are volume fractions of Ω locally occupied by the biofilm, u, and its

172

components- we denote by H, V, C and Y the components Hmp, N orV ,

173

Cyd and complex Cyd−N O• respectively. The dissolved substrates are oxy-

174

gen, glucose, nitric oxide and DPTA NONOate which are described in terms

175

of their concentrations N1 [µM ], N2 [µM ], N3 [µM ] and N4 [µM ] respectively. 8

176

Using the above notation, we obtain the following model:

                                                                          

∂u ∂t

N2 N1 ) u− = ∇(D(M )∇u) + µu (1.82 + κOu + N1 κG + N2 {z } | growth of biofilm

∂H ∂t

N3 β u κ+N | {z 3 }

(1a)

inactivation of biofilm

N1 N3 u − dH H = ∇(D(M )∇H) + µH | {z } κOH + N1 κN H + N3 | {z } natural decay

(1b)

activation of Hmp

∂V ∂t

κ2 N3 = ∇(D(M )∇V ) + µV 2 OV 2 u − dV V |{z} κOV + N1 κN V + N3 | {z } natural decay

(1c)

activation of NorV

∂C ∂t

= ∇(D(M )∇C) −

N4 µC 4 3 4 C κ +N | N C{z 3 }

+

transformation to complex

∂Y ∂t

= ∇(D(M )∇Y ) +

N4 µC 4 3 4 C κ +N | N C{z 3 }

(1d)

back transformation from complex

activation after binding N O• to C

9

κ4 dC 4 N Y 4 Y κ +N | N Y{z 3 }

κ4 − dC 4 N Y 4 Y κ +N | N Y{z 3 }

back transformation to C

(1e)

                                                                                                                           177

∂N1 ∂t

N1 N3 N1 1−γ H − ν1 C = ∇(dN1 ∇N1 ) − µH M0 γ κN H + N3 κOH + N1 κr + N1 {z } | {z } | respiration

uptake by Hmp

1 N1 N2 N1 1−γ u− α N3 − µu M0 γ κOu + N1 κG + N2 2 κo + N1 {z } | {z } | Autoxidation

uptake by biofilm

(1f )

∂N2 ∂t

= ∇(dN2 ∇N2 ) −

N2 µu M0 N1 ) u (1.82 + γ κOu + N1 κG + N2 | {z }

(1g)

uptake by biofilm

∂N3 ∂t

= ∇(dN3 ∇N3 ) −

N3 N1 µH M0 H γ κN H + N3 κOH + N1 | {z } degredation by Hmp

−

µV M0 N3 κ2OV V γ κN V + N3 κ2OV + N12 {z } | degradation by NorV

N4 − µC M0 4 3 4 C + κN C + N3 | {z } binding to C

κ4 dC M0 4 N Y 4 Y κN Y + N3 | {z }

back transformation from complex

N1 N3 + µN N4 −α | {z } κo + N1 | {z } production of nitric oxide

(1h)

Autoxidation

∂N4 ∂t

= ∇(dN4 ∇N4 ) −

ν2 N4 | {z }

DPTA dissociation

(1i) (1)

Here M := u + H + V + C + Y is the total volume occupied by the 10

178

bacterial biomass. The diffusion of biomass is described by a nonlinear den-

179

M sity dependent coefficient D(M ) = δ (1−M [m2 d−1 ] in which δ [m2 d−1 ] is a )4

180

positive parameter and smaller than the diffusion coefficient of dissolved sub-

181

strates in the fluid [7]. The diffusion coefficient D(M ) displays two nonlinear

182

effects: (i) a porous medium degeneracy, i.e. D(M ) vanishes as M ≈ 0 and

183

(ii) a fast diffusion singularity as M approaches unity. The porous medium

184

degeneracy, M 4 , guarantees that the biofilm/water interface propagates with

185

finite speed and does not spread significantly if the biomass density is small,

186

0 < M  1. It also results in the formation of a sharp interface between

187

biofilm and surrounding liquid (i.e. initial data with compact support lead

188

to solutions with compact support). The second effect (ii) at 0  M < 1

189

enforces the condition that the solution is bounded by unity [7, 48]. This is

190

counteracted by the degeneracy as M = 0 at the interface. Consequently,

191

M squeezes in the biofilm region and approaches its maximum value of 1.

192

Hence, the interaction of both nonlinear diffusion effects with the growth

193

term is necessary to describe spatial biomass spreading [7]. To make the

194

interpretation of the the relevant interactions introduced in the model easier,

195

a rough sketch of the interactions between species and substrates is shown

196

in Figure 1.

4

197

198

Below is a description of each reaction term included in each equation of system (1):

199

• (1a). Since the growth of the biofilm u depends on the available glucose

200

and is enhanced in the aerobic regime (c.f. [37]), it is described by 11

Figure 1: Sketch of the interactions between bacteria (in biofilms), oxygen, N O•, nutrient and N O• donor; solid lines denote a positive influence(production), dashed lines a negative influence(consumption)

201

standard and dual Monod kinetic terms in which κOu [µM ] and κG [µM ]

202

are the oxygen and glucose half saturation concentrations respectively

203

and µu [h−1 ] is the maximum growth rate. The constant number 1.82

204

and values of parameters κG [µM ] and µu [h−1 ] are found by fitting the

205

experimental data provided in the supplementary material of [41] with

206

the Monod function. Eradication of the biofilm due to exposure to

207

N O• is captured by a Monod kinetic term with maximum decay rate

208

β[h−1 ] and half saturation concentration κ[µM ].

209

• (1b), (1c). Species H is produced in the presence of N O•[µM ] and oxy-

210

gen while the production of component V is stimulated by N O • [µM ]

211

under anaerobic contribution. This leads to the introduction of the 12

212

standard dual Monod kinetic term for H with half saturation concen-

213

trations κN H , κOH and maximum production rate µH . While produc-

214

tion of V is described by an inhibition Hill function with exponent

215

2, half saturation concentrations κN V , κOV and maximum production

216

rate µV . The inhibition term is described by a hill function to give

217

a more threshold-like behavior to the production of V . It is assumed

218

that H and V decay naturally and (for simplicity) have the same lysis

219

rate.

220

• (1d), (1e). An important impact of N O• on bacterial cells is its in-

221

hibition effect on respiration activity by reversibly binding Cyd [28].

222

Under normal respiratory condition, oxygen binds a heme at the active

223

site of E.coli’s Cyd and is released as H2 O following a two-electron

224

reduction [33]. As the affinity of Cyd for N O• is much greater than

225

for oxygen, even low N O• concentration can inhibit the activity of

226

Cyd [29]. In the proposed model, reversible binding of N O• to Cyd

227

which acts as a temporary N O• sink and indicator of respiratory in-

228

hibition is described by a Hill function with exponent 4. The reason

229

for choosing exponent 4 is the sensitivity of Cyd activity to N O• that

230

needs even more threshold-like behavior. Moreover it is shown in [31]

231

that to reflect the biologically observed oscillatory phenomena in pro-

232

tein interaction, switch functions should be Hill functions with expo-

233

nent n ∈ {2, ..., 19}. It is proved mathematically in [31] that by these

234

exponents and by making the transcription time small enough, one is 13

235

always able to obtain damped oscillatory behaviour. Half saturation

236

concentrations, production and back transformation rates are shown

237

by κN C , κN Y and µC , dC respectively.

238

• (1f). Oxygen consumption by Hmp and u is described by dual Monod

239

kinetics at a rate proportional to the production and growth rates of

240

H and u, yield coefficient γ[−] and maximum cell density M0 [gm−3 ].

241

Respiration activity of Cyd which degrades the oxygen concentration

242

and occurs in the absence of nitric oxide is described by a Monod func-

243

tion with half saturation concentration κr and uptake rate ν1 . The

244

model also includes the interaction between N O• and oxygen i.e au-

245

toxidation. Many damaging effects of N O• result from the activity of

246

its autoxidation products. In our model, autoxidation which degrades

247

oxygen and nitric oxide is described by a Monod kinetic term with

248

κo [µM ] as the half saturation concentration. Considering the reaction

249

2N O • +O2 → 2N O2 •, degradation rates for N O • [µM ] and oxygen

250

are α[h−1 ] and 21 α[h−1 ] respectively.

251

• (1g). Nutrient consumption by the biofilm is described by standard and

252

dual Monod kinetics at a rate proportional to the growth rate µu [h−1 ],

253

yield coefficient γ[−] and maximum cell density M0 [gm−3 ].

254

• (1h). N O• degradation by H, V and C is described by the same kinet-

255

ics that were described before at rates proportional to the production

256

rate, maximum cell density and yield coefficient. N O • [µM ] is gener14

257

ated in the fluid after adding DPTA. The maximum rate for N O •[µM ]

258

production, µN [h−1 ], is defined as µN =

2Ln(2) −1 [h ] τ1 2

in which τ 1 [h] is 2

the half-life time of DPTA respectively.

259

• (1i). DPTA is degraded due to dissociation. Its degradation rate is

260

calculated as ν2 =

261

Ln(2) −1 [h ]. τ1 2

The parameters and their values used in computer simulation are sum-

262

263

marized in Table 1.

264

Remark 2.1. Since no appropriate experimental data exists to determine

265

the values of parameters for the proposed model, we have adapted the biofilm

266

parameter values from [6, 13, 15] and have tried to estimate N O• related

267

reactions based on realistic assumptions from [37–41] for well-mixed cultures.

268

3

Simulation set-up

269

3.1

Boundary condition

270

To study the mathematical model (1) and not be hindered by numerical

271

artifacts, we restrict ourselves to a simple square domain Ω = [0, L] × [0, L]

272

and place one semi-spherical colony at the center of the substratum y = 0.

273

We assume that the substratum is impermeable to biomass and dissolved

274

substrates. At the lateral boundaries, x = 0 and x = L, a homogeneous

275

Neumann condition is applied for the dependent variables, assuming that

276

the domain is part of a continuously repeating larger domain. We assume

277

that oxygen and glucose enter the system through the top boundary (i.e. 15

Table 1: Model parameters for system (1) used for computer simulation Parameter Symbol Growth rate of u µu Production rate of H µH Production rate of V µV Transformation rate from C to Y µC Back transformation rate from Y to C dC Oxygen monod half saturation for u κOu Nutrient monod half saturation for u κG N O• monod half saturation for H and V κN H,N V Oxygen monod half saturation for H κOH Oxygen monod half saturation for V κOV N O• monod half saturation for C κN C Oxygen monod half saturation for κN Y complex formation Oxygen monod half saturation κr for respiration Oxygen monod half saturation κo for autoxidation Decay rate of H and V dH,V Decay rate of u by N O• β Respiration rate ν1 N O• autoxidation rate α Maximum cell density M0 Oxygen diffusion coefficient dN1 Nutrient diffusion coefficient dN2 N O• diffusion coefficient dN3 DPTA diffusion coefficient dN4 Concentration boundary layer thickness λ yield coefficient γ Nutrient bulk concentration (N2 )∞ Half-life time of DPTA τ1 2

16

Value 0.07 1 0.5 1000 1000 3.986 2.6 0.01 25 0.1 2 0.1

Dimension h−1 h−1 h−1 h−1 h−1 µM µM µM µM µM µM µM

2

µM

100

µM

0.25 0.01 117 × 105 15 6500 8.3 × 10−6 0.04 × 10−14 12 × 10−6 4.1 × 10−6 0.5 0.63 1170 2.5

h−1 h−1 µM h−1 h−1 µM m2 h−1 m2 h−1 m2 h−1 m2 h−1 mm [−] µM h

278

y = L) so that a Robin condition is imposed for oxygen and glucose over this

279

segment. It is also assumed that DPTA NONOate is initially available in the

280

domain and releases N O • [µM ] and that the top boundary is impermeable

281

to them. Hence, a homogeneous Neumann condition is imposed for both

282

N O • [µM ] and DPTA at y = L. Thus, the boundary condition for each

283

species and substrate is as follows:  ∂n u = ∂n H = ∂n V = ∂n C = ∂n Y = 0, at ∂Ω         N1 + λ∂n N1 = (N1 )∞ , at y = L      N2 + λ∂n N2 = (N2 )∞ , at y = L       ∂n N2 = ∂n N1 = 0, at x = 0, L and y = 0        ∂n N3 = ∂n N4 = 0, at ∂Ω

(2)

284

where (N1 )∞ and (N2 )∞ are the bulk concentrations of oxygen and glu-

285

cose respectively. Parameter λ[mm] can be interpreted as the concentration

286

boundary layer thickness and ∂n denotes the outward normal derivative. The

287

concentration boundary layer mimics the contribution of the convective con-

288

tribution of external bulk flow to substrate supply. A small bulk flow veloc-

289

ity implies a thick concentration boundary layer, while a thin concentration

290

boundary layer represents fast bulk flow [8].

291

3.2

292

The two dimensional computational domain Ω is divided into two regions

293

Ω1 (t) = {(x, y) ∈ Ω ⊂ R2 : M (x, y, t) = 0} that describes the aqueous phase

Initial condition

17

294

(bulk liquid, channels and pores of a biofilm) without biomass, and region

295

Ω2 (t) = {(x, y) ∈ Ω ⊂ R2 : M (x, y, t) > 0}, which is the actual biofilm with

296

positive biomass density. Both regions are separated by the biofilm/water

297

¯ 1 (t) ∩ ∂ Ω ¯ 2 (t). The biofilm region Ω2 (t) usually consists of interface Γ(t) = ∂ Ω

298

several colonies which might merge as they expand. In fact, Ω1 and Ω2 and

299

their interface are time dependent because of changing the biofilm structure

300

due to growth.

301

One semi-spherical colony is placed at the center of the substratum. The

302

initial values of biomass volume fractions u and C are u = 0.9 and C = 10−4

303

and the initial values of the other biomass volume fractions are set to zero.

304

The values of oxygen and glucose are set initially to their bulk concentrations

305

and the default value of DPTA at t = 0 is set to 50[µM ] according to the

306

experimental data presented in [41]. DPTA is present initially in the liquid

307

phase, releases N O• in this region and both DPTA and N O• are transported

308

into the biofilm by diffusion.

309

3.3

310

Each of the two nonlinear diffusion effects (i) and (ii) raise unique chal-

311

lenges for numerical simulation, see [13,15,16] for more details. It was shown

312

in [13, 14, 16] that error-controlled adaptive time integration methods can

313

overcome the difficulties that methods with fixed time-steps experience in

314

solving singular ODE systems (obtained by spatial discretization of biofilm

315

models with nonlinear effects (i) and (ii), by standard Finite Volume meth-

Numerical Method

18

316

ods). Here also we use a finite volume scheme to discretize our PDE model

317

and apply the third order embedded Rosenbrock-Wanner method with 4

318

stages (ROS3PRL) [34] (the same time integration method as in [13,14,16]),

319

to solve the semi-discrete problem. For a more detailed description of the

320

numerical method we refer to [13].

321

For the numerical treatment we non-dimensionalize model (1): x˜ = x/L

322

and t˜ = t × µH for the independent variables, where L is the height of the

323

computational domain and scaling factor for time is µH [h−1 ]. Choosing ei-

324

ther µu = 0.07[h−1 ] or µC = dC = 1000[h−1 ] as the time scaling factor

325

results in a computational error due to the choice of characteristic time scale

326

(i.e. very small or large values respectively). To avoid this error, we chose

327

µH as the time scaling factor. The substrate concentration variable N2 is

328

˜2 = non-dimensionalized with respect to its bulk concentration N

329

˜4 = DPTA is non-dimensionalized with respect to its initial value N

330

The scaling factors to non-dimensionalize N1 (O2 ) and N3 (N O•) are cho-

331

sen as 10[µM ] and 1[µM ] for simplicity. Although we use the dimensionless

332

formulation for the numerical treatment, we will make our choices of param-

333

eters based on the dimensional values, and discuss our subsequent simulation

334

experiments in terms of the dimensional values, for ease of physical interpre-

335

tation.

19

N2 (N2 )∞

and

N4 . N4 (0)

336

3.4

Postprocessing and quantities of interest

For better interpretation of the results of computer simulations, we display the time course of N O• concentration and the spatial distribution of N O• in the computational domain through two-dimensional visualizations. To quantify the impact of changing parameters on N O• detoxification, we define the following lumped output variables as quantities of interest:

 R I(t) = I(x, y, t) dx dy, I = u, H, V, N3  Ω2      R R  N1  3  N O • detox. by H = t Ω2 µHγM0 κN HN+N H(x, y, t) dx dy dt,  3 κOH +N1      R R µV M0 N3  κ2OV     N O • detox. by V = t Ω2 γ κN V +N3 κ2OV +N12 V (x, y, t) dx dy dt,

R R  κ4 Y   N O • degradation by C = t Ω2 dC M0 κ4 N+N 4 Y (x, y, t)   3 NY       N34  −µC M0 κ4 +N  4 C(x, y, t) dx dy dt,  3 NC       Autoxidation = R R α N1 N (x, y, t) dx dy dt t Ω κo +N1 3 337

Remark 3.1. Besides the above degradation pathways, N O• concentration

338

decreases in the domain due to evaporation i.e. through the gas phase. It is

339

shown in [37, 41] that the amount of N O• which leaves the system as a gas

340

is not negligible. Despite this fact, since we do not consider the contribution

341

of convective flux in transporting the dissolved substrates, we do not have to

342

take this pathway into account in our current study. 20

343

4

Results

344

4.1

345

4.1.1

An Illustrative Simulation Temporal behavior of N O• (A)

(B)

1.6

(C)

16

1.4

7

14

(N1)∞=0 [µΜ]

6

(N1)∞=5 [µΜ]

1.2

12

1

10

(N1)∞=10 [µΜ]

5

0.8

No

No

No

4 8

3 0.6

6

0.4

4

0.2

2

0

0

1

2

3

4

0

5

2

1

0

1

2

t

3

4

0

5

(D)

1

2

3

5

(F)

2.5

(N1)∞=15 [µΜ]

4

t

(E)

3.5

3

0

t

0.7

0.6

(N1)∞=20 [µΜ]

2

2.5

(N1)∞=50 [µΜ]

0.5 1.5

1.5

No

0.4

No

No

2

0.3

1

1

0.2 0.5

0.5

0

0.1

0

1

2

3

4

5

0

0

1

2

t

3

t

4

5

0

0

1

2

3

4

t

Figure 2: Time course of N O • (N3 (t)) concentration in the biofilm region for different values of oxygen bulk concentration.

346

Figure 2 illustrates the interplay between oxygen concentration, activity

347

of N O• reductants and N O• concentration within the biofilm under three

348

conditions: aerobic, microaerobic and anaerobic regimes. As shown in Fig-

349

ures 2A-F, the concentration of N O• initially increases since there is no Hmp

350

and N orV at the beginning to reduce N O•. Reduction of N O• through au-

351

toxidation is not significant. Despite the increase in the N O• concentration

352

at the initial stage, it is not large enough to inhibit respiration, indicating

353

slow reduction in oxygen concentration at this step. By increasing N O• con21

5

354

centration and exceeding 0.5[µM ], respiration is inhibited and oxygen levels

355

go up. Up until a specific time which depends on the oxygen bulk concentra-

356

tion, either Hmp or N orV are produced gradually. After that, N O• detoxi-

357

fication occurs which moves its concentration down. At this phase, there is

358

not enough N O• to inhibit respiration. So, oxygen outcompetes N O• and

359

respiration activity resumes which consequently decreases the oxygen concen-

360

tration and limits all oxygen dependent activities of biomass species. Note-

361

worthy here is the oscillatory behavior of N O• in the microaerobic regime in

362

Figures 2B-E. The reason for this phenomenon which prolongs the time for

363

N O• clearance, as observed experimentally for the well mixed culture [41],

364

is due to interactions between oxygen dependent activity of Hmp and Cyd

365

and N O• consumption by Hmp. First, N O• concentration decreases within

366

the biofilm due to consumption by Hmp. Consequently, respiration activity

367

resumes which reduces the oxygen concentration and deactivates Hmp. This

368

results in an increase in the N O• concentration levels which subsequently

369

inhibits the respiration activity. The oscillatory behavior lasts for approx-

370

imately 30 min and then N O• concentration is totally diminished in the

371

biofilm. The same behavior (as in the anaerobic regime) is observed in Fig-

372

ure 2F for (N1 )∞ = 50[µM ], however N O• reaches to submicromolar regime

373

a bit faster.

374

To further assess the predictive power of the model, we investigate the

375

role of Hmp in removing N O• and effect of gene mutation on the bacterial

376

defense system under an aerobic condition. For this purpose, we simulate 22

(A)

(B) 1.0 0.8 0.6

40

0.4

20

0.2

0 0

1

2

t

3

4

5

0.0

5 Autoxidation Hmp NO

80 NO consumption

60

100

NO

NO consumption

80

Autoxidation Hmp NO

60

3

40

2

20

1

0

0

1

2

t

3

4

Figure 3: Distribution of N O• consumption via autoxidation and detoxification by Hmp for the culture (left) with Hmp and (right) without Hmp. The bulk concentration of oxygen is set to (N1 )∞ = 50 [µM ].

377

the distribution of N O• consumption in a system with and without Hmp.

378

The bulk concentration of oxygen is set to (N1 )∞ = 50 [µM ]. The results

379

which are displayed in Figure 3 highlight the effects of Hmp in removing

380

N O• and consequently its deficiency in bacterial eradication. Although the

381

N O• curve for both the mutant (without Hmp) and wildtype (with Hmp)

382

cultures are initially very close to each other, at t ≈ 0.2[h] they start to

383

diverge due to the production of Hmp and its dominant impact on the N O•

384

detoxification in the aerobic regime. As shown in Figure 3, the model predicts

385

that in a wild-type culture enzymatic N O• detoxification by Hmp results in

386

a quick degradation of N O• while in a mutant culture N O• concentration

387

gradually decreases via autoxidation, hence it can remove bacterial biofilm

388

efficiently. The obtained computational results are in good agreement with

389

the experimental data provided in [37] 23

4

NO

100

0 5

390

4.1.2

Spatial distribution of N O•

t = 0.3[h]

t = 0.68[h]

t = 0.81[h]

t = 0.98[h]

Figure 4: Snapshot of N O • (N3 ) and oxygen (N1 ) concentration at t = 0.3, 0.68, 0.81, 0.98. The color coding refers to the N O• concentration and greyscale isolines show the oxygen concentration. The biofilm interface is indicated in red. Oxygen bulk concentration is N1 = 15[µM ].

391

To explore the spatial distribution and variation of N O• in the computa-

392

tional domain, we provide a two dimensional visualisation and compute the

393

average of N O• and Hmp along with standard deviation in Figures 4-6. The

394

latter quantities of interest are computed as: 24

t = 1.06[h]

t = 1.22[h]

t = 1.33[h]

t = 3[h]

Figure 5: Snapshot of N O• and oxygen (N1 ) concentration at t = 1.06, 1.22, 1.33, 3. The color coding refers to the N O• concentration and greyscal isolines show the oxygen concentration. The biofilm interface is indicated in red. Oxygen bulk concentration is (N1 )∞ = 15[µM ].

25

    

AveN3 (y) := SDN3 (y) :=

1 L

RL 0

N3 (x, y, t∗ )dx

q R 1 L L

0

(3) (N3

(x, y, t∗ )

2

− AveN3 (y)) dx

395

Average and standard deviation of Hmp are computed in the same way as

396

in (3).

397

The spatial distribution of N O• at different times under a microaerobic

398

regime is shown in Figures 4-5. The color coding refers to the N O• con-

399

centration and greyscale isolines show the oxygen concentration. Graphical

400

results presented in Figures 4-5 reveal that production of N O• is initiated

401

from the liquid phase near the top boundary nevertheless the gradient of its

402

concentration between the biofilm and aqueous phase is not initially steep

403

due to fast diffusion of N O•. At some time greater than 0.3, enough Hmp is

404

produced to remove N O•. Thus at t = 0.68, the N O• concentration within

405

the biofilm is less than 0.3 while it is larger than 1 in the liquid phase. At

406

this time, there is a small amount of N O• to bind Cyd so respiration occurs

407

and oxygen concentration decreases which disables Hmp and consequently

408

prevents N O• detoxification. This leads to a slight increase in the N O• con-

409

centration at t = 0.81, c.f. Figure 4. This biologically oscillatory behavior in

410

the N O• concentration continues at t = 1.06, 1.22 and t = 1.33 and finally

411

stops at t = 3. At this time N O• concentration reaches to its minimum

412

value everywhere in the domain and later on is diminished, see Figures 4-5.

413

An important finding is that for the considered geometry, the observed os-

414

cillations are uniform in space (i.e. this is a temporal phenomenon and does 26

415

not give rise to a spatial pattern). 0.8

0.14

0.7

0.12 0.1

0.6

0.08

H(y)

NO(y)

0.5

0.4

0.06 0.04

0.3 0.02 0.2

0

0.1

0

-0.02

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

-0.04

1

y

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

y

Figure 6: The average of N O • (N3 (y)) and biomass volume fraction Hmp at t = 0.5[h]. Oxygen bulk concentration is (N1 )∞ = 50[µM ]. 416

The average of N O• and Hmp along with the standard deviation are

417

plotted in Figure 6 to show their spatial deviation. The results are obtained

418

under aerobic conditions at t = 0.5[h]. We observe that maximum varia-

419

tions in the mean value of N O• and Hmp occur near the substratum and

420

become smaller towards the top boundary. The results in Figure 6 indi-

421

cates the heterogeneous distribution of N O• and Hmp within the biofilm

422

and homogeneous distribution of N O• in the liquid phase.

423

4.2

424

Biofilm defense vs initial value and half-life time of N O• donor

425

An underlying assumption of our work is that N O• is generated by DPTA

426

NONOate exogenously. To examine the effect of initial concentration of N O•

427

donor and its release rate on the N O• detoxification, we consider two scenar-

428

ios. In the first case it is assumed that one type of N O• donor with different 27

(A) 35

30

(N4)0=50 [µM] (N4)0=100 [µM] (N4)0=200 [µM] (N4)0=300 [µM] (N4)0=500 [µM]

25

NO

20

15

10

5

0

0

1

2

3

4

5

t

(B)

(C) 100%

90%

90%

%NO consumption via Autoxidation

100%

% NO consumption by Hmp

80% 70%

(N4)0=50 [µM] (N4)0=100 [µM] (N4)0=200 [µM] (N4)0=300 [µM] (N4)0=500 [µM]

60% 50% 40% 30% 20% 10% 0%

(N4)0=50 [µM] (N4)0=100 [µM] (N4)0=200 [µM] (N4)0=300 [µM] (N4)0=500 [µM]

80% 70% 60% 50% 40% 30% 20% 10%

0

1

2

3

4

0%

5

0

1

2

t

3

4

5

t

Figure 7: (A): Time course of N O• (N3 (t)) concentration for different initial values of DPTA. (B),(C): N O• consumption by Hmp and autoxidation for different initial values of DPTA. The results are for (N1 )∞ = 50[µM ]. 10

100%

9

90%

t1/2=0.5 [h] t1/2=2.5 [h]

8

80%

% NO consumption

7

NO

6 5 4 3 2

60%

Consumption by Hmp, t1/2=0.5[h] Autoxidation, t1/2=0.5[h] Consumption by Hmp, t1/2=2.5[h] Auoxidation, t1/2=2.5[h]

50% 40% 30% 20%

1 0

70%

10%

0

1

2

3

4

0%

5

t

0

1

2

3

4

5

t

Figure 8: Time course of N O • (N3 (t)) concentration and the percentage of its consumption by Hmp and autoxidation for large and small value of donor release rate. The results are for (N1 )∞ = 50[µM ] and (N4 )0 = 50[µM ]. 28

429

initial concentrations is used to generate N O•. While in the second scenario,

430

two different types of DPTA with different half-life times and same initial

431

concentrations are added to the system. As the results in Figure 7 show,

432

increasing the initial concentration of N O• donor provides more N O• to the

433

system which can not be detoxified by the available reducing agents and con-

434

sequently its clearance time increases. However, limitless increasing of DPTA

435

is not practically possible as it might be toxic. The distribution of the N O•

436

consumption pathways indicates that by increasing the initial concentration

437

of N O• donor, Hmp becomes the minor reductant and N O• is degraded

438

mostly via autoxidation, see Figure 7B,C. For higher initial concentrations

439

of DPTA, interaction between N O• and oxygen increases in the liquid re-

440

gion which results in N O• neutralization through autoxidation. Hence, only

441

a small fraction of the generated N O• reaches the biofilm that is detoxified

442

by stimulated enzymes. In the second scenario, two different types of N O•

443

donors with different half-lives, are administered. By increasing the release

444

rate i.e. decreasing the half-life time, more N O• is released which can not be

445

detoxified fast by defense agents and takes longer to be cleared, (c.f Figure

446

8). As more N O• is generated in the liquid phase (by increasing the release

447

rate of the N O• donor), the activity of autoxidation increases which results

448

in the substrate limitation for Hmp. This consequently reduces the ability of

449

Hmp to participate in N O• detoxification during the first hour, see Figure 8.

450

When Hmp activity is restored after ∼ 30 minutes, most of the N O• has al-

451

ready been consumed through autoxidation. As a result, it is predicted that 29

452

enzymatic detoxification of N O• in the wild-type culture will be decreased.

453

This finding indicates that there is not a great difference between clearance

454

time of N O• in the wild-type and mutant culture (data not shown).

455

In summary, generating N O• with large amounts of N O• donor with

456

short half-life is more efficient in the sense that more N O• is released, the

457

activity of Hmp is reduced and the time period required for N O• to ap-

458

proach the submicromolar regime is increased. All these can increase the

459

inhibitory effect of N O• as an antimicrobial agent on biofilm growth. How-

460

ever, we should note that activity of N O• autoxidation products including

461

nitrogen dioxide (N O2 •), nitrous anhydride (N2 O3 ) and dinitrogen tetroxide

462

(N2 O4 ) results in many damaging effects indicating that autoxidation should

463

be brought under control [24].

464

4.3

465

Our objective in this section is to assess the effect of the initial size of a

466

biofilm colony, that might change the N O• flux and produced Hmp, on the

467

N O• detoxification. For this purpose, we change the initial radius of the

468

biofilm colony from r = 0.2 [mm] to r = 0.04 [mm] and compute the total

469

amount of N O• and produced Hmp within the biofilm region in the aerobic

470

and microaerobic regime, results are shown in Figure 9. We observe that by

471

decreasing the initial size of the biofilm colony, concentration of N O• inside

472

the biofilm increases under both aerobic and microaerobic conditions which

473

consequently increases its clearance time. By decreasing the initial size of

N O• detoxification vs initial size of biofilm colony

30

(A)

(B)

4

14

Small colony Large colony

3.5

Small colony Large colony

12

3

10

2.5

NO

NO

8 2

6 1.5 4

1

2

0.5

0

0

2

4

6

8

0

10

0

2

4

6

t 0.02

0.02

0.018

0.018

0.016

0.014

Large colony Small colony

0.012

0.01

H

H

0.012

0.01

0.008

0.008

0.006

0.006

0.004

0.004

0.002

0.002

0

2

10

0.016

Large colony Small colony

0.014

0

8

t

4

6

8

0

10

0

2

4

t

6

8

10

t

0.4

0.3

Small colony Large colony

0.35

Small colony Large colony

0.25 0.3 0.2

0.2

u

u

0.25

0.15

0.15 0.1 0.1 0.05 0.05

0

0

2

4

6

8

0

10

t

0

2

4

6

8

10

t

Figure 9: Time course of N O • (N3 (t)) concentration, produced Hmp and biofilm concentration u(t) for (A) aerobic, (N1 )∞ = 50[µM ], and (B) microaerobic, (N1 )∞ = 10[µM ], regime. The initial radius of biofilm colony is chosen as r = 0.2 [mm] for large biofilm and r = 0.04 [mm] for small biofilm. The results are for (N4 )0 = 50[µM ].

31

474

the biofilm, less N O• diffuses into the biofilm and due to the production of

475

less Hmp its enzymatic detoxification also decreases. However, reduction in

476

the amount of incoming N O• is less than that in its detoxification by Hmp,

477

indicating the existence of more N O• inside the biofilm with smaller initial

478

size. As the results in the third row of Figure 9 show, N O• can eradicate

479

the biofilm in a more effective way if biofilm size is small. Noteworthy here

480

is the absence of oscillatory behavior of N O• in the microaerobic regime,

481

(N1 )∞ = 10[µM ]. By decreasing the initial size of biofilm colony, less Hmp

482

is produced, c.f the second row of Figure 9, which can not compete with Cyd

483

for oxygen consumption and this prevents oscillation.

484

In summary, the defense system of a thin biofilm against N O• is weaker

485

compared with that of a thick biofilm. This increases N O• level which is

486

able to eradicate the biofilm. The reason could be either reduction in the

487

amount of N O• that diffuses into the biofilm and its consumption by Hmp

488

or production of lower amounts of Hmp or a combination of both. The

489

simulation results suggest that treatment of bacterial stain by N O•, at the

490

early stage of its development, is more efficient in the sense of inhibiting the

491

virulence factors.

492

4.4

493

The effect of nutrient deprivation on N O• deactivation

494

Pathogens use reducing agents, whose activity is stimulated by N O• as a

495

stressor, in order to ward-off N O•. The presence of additional stressors such 32

(A)

(B)

0.7

4

(N2)∞=1170 [µM] (N2)∞=1.17 [µM]

(N2)∞=1170 [µM] (N2)∞=1.17 [µM]

3.5

0.6

3

0.5

2.5

NO

NO

0.4 2

0.3 1.5 0.2

1

0.1

0

0.5

0

1

2

3

4

0

5

0

1

2

t

3

4

5

t

(C)

(D)

3.5

10

(N2)∞=1170 [µM] (N2)∞=1.17 [µM]

(N2)∞=1170 [µM] (N2)∞=1.17 [µM]

9

3 8 2.5

7 6

NO

NO

2

1.5

5 4 3

1

2 0.5 1 0

0

1

2

3

4

0

5

t

0

1

2

3

4

5

t

Figure 10: Time course of N O • (N3 (t)) concentration in the aerobic, (N1 )∞ = 50[µM ], and microaerobic, (N1 )∞ = 15[µM ], regime for limited and unlimited nutrient concentration. Panel (A) shows the result for thick biofilm with initial radius r = 0.2 [mm] under the aerobic condition. Panel (B) is for thin biofilm with initial radius r = 0.04 [mm] in the aerobic regime. Panel (C) reports the result for thick biofilm in the microaerobic regime and panel (D) shows the result for thin biofilm under the microaerobic condition. The results are for (N4 )0 = 50[µM ].

33

(A)

(B)

1.6

(N2)∞=1170 [µM] (N2)∞=1.17 [µM]

1.4

1.4 1.2

(N2)∞=1170 [µM] (N2)∞=1.17 [µM]

1.2 1

N1

N1

1

0.8

0.6

0.6

0.4

0.4

0.2

0.8

0.2

0

1

2

3

4

0

5

t

0

1

2

3

4

5

t

Figure 11: Time course of oxygen (N1 (t)) concentration in the microaerobic, (N1 )∞ = 15[µM ], regime for limited and unlimited nutrient concentration. Panel (A) is for thick biofilm with initial radius r = 0.2 [mm] and panel (B) is for thin biofilm with initial radius r = 0.04 [mm]. The results are for (N4 )0 = 50[µM ]. 496

as nutrient deprivation and acid stress can affect the mechanism of action

497

leading to detoxification [17]. Biofilm growth depends upon the availability

498

of carbon sources, and increases under aerobic conditions. Nutrient depri-

499

vation limits the production of u which reduces the formation of Hmp. On

500

the other hand, it leads to an increase in the oxygen concentration which

501

increases the value of the biomass volume fraction Hmp. This nonlinear

502

interaction between nutrient supply and N O• detoxification makes it com-

503

plicated to predict how starvation affects the biofilm defense system against

504

N O•. The aim in this section is to explore the influence of nutrient limita-

505

tion on N O• deactivation under aerobic and microaerobic regimes for thin

506

and thick biofilms. For this purpose we change the glucose bulk concentra-

507

tion and compute the total amount of N O• in the biofilm under limited and

508

unlimited nutrient conditions. 34

509

The results shown in Figures 10A,B reveal that under the aerobic con-

510

ditions, lack of nutrient does not considerably change the amount of N O•

511

and its clearance time within the biofilm. Although decreasing the nutrient

512

supply, decreases biofilm growth nevertheless, over the clearance time inter-

513

val, this term is not dominated by N O• dependent decay term of biofilm.

514

Therefore, biofilm volume fraction, u, does not decrease significantly thus it

515

does not have any impact on the activity of Hmp. This highlights the role

516

of oxygen on the production of Hmp for the considered simulation set-up.

517

As oxygen is not limiting in the aerobic and even microaerobic regimes, tem-

518

poral behavior of N O• under limited and unlimited nutrient conditions is

519

the same. Results in panels 10C,D show the N O• concentration in a thick

520

and thin biofilm respectively, under microaerobic conditions for limited and

521

unlimited nutrient supply. We observe that, compared to the aerobic regime

522

(Figures 10A) the clearance time of N O• in the thick biofilm increases in the

523

microaerobic regime regardless of the amount of carbon source. This pro-

524

vides enough time for the N O• mediated decay term of biofilm to dominate

525

its growth term. Under the nutrient deprivation conditions, this leads to an

526

increase in the oxygen concentration (see panel (A) of Figure 11) which sub-

527

sequently increases the production of Hmp and lowers the maximum value

528

of N O•, c.f Figure 10C. However, the trend in the thin biofilm is completely

529

different. When the biofilm size is small, the difference between oxygen con-

530

centration under limited and unlimited nutrient conditions is not significant

531

(c.f panel (B) of Figure 11), so production of Hmp is mainly controlled by u. 35

532

Clearance time of N O• within a thin biofilm under microaerobic conditions

533

is large enough to observe the effect of decay of u (due to starvation) on

534

reducing the production of Hmp and making the N O• detoxification slow.

535

Our results suggest that depending on the biofilm size and oxygen con-

536

centration levels, starvation can change the defense system against N O•.

537

Under aerobic conditions, the impact of starvation on the removal of N O•

538

is not significant. Nevertheless, thin biofilms have weaker defense systems

539

under nutrient limitation stress although deviation is very small. Under a mi-

540

croaerobic regime and upon applying nutrient deprivation stress, it is better

541

to treat the biofilm before it gets thick otherwise the defense system becomes

542

very strong and N O• is deactivated fast.

543

The effect of nutrient deprivation stress on a pathogen defense system is

544

studied experimentally in [17]. It is shown in [17] that under limited nutrient

545

conditions, N O• detoxification by Hmp fails in aerobic E. coli cultures. The

546

results that we obtained for thin biofilms in the microaerobic regime are in

547

good agreement with the results presented in [17] for replete carbon source.

548

The reasons that we got the same results but under different condition could

549

be simulation set-up, boundary conditions and parameter values in the re-

550

action kinetics. In order to make the simulation set-up consistent with the

551

experimental condition in [17] and show the effect of deprivation stress on

552

the failure of N O• detoxification by Hmp more clearly, we repeated the sim-

553

ulations for a thin biofilm in the microaerobic regime with a Robin boundary

554

condition for DPTA. In the new simulation set-up, DPTA enters to the sys36

555

tem continuously, while there is not enough Hmp to defend against N O•,

556

see Figures 12A,B. As the results show, under starvation stress, the viru-

557

lence factor (detoxification by Hmp) is inhibited and more N O• exists in

558

the system which can effectively remove bacterial biofilm, c.f. Figure 12C. A

B 0.012

(N2)∞=1170[µM] (N2)∞=1.170 [µM]

10

0.01

8

0.008

6

0.006

H

NO

12

4

0.004

2

0.002

0

0

2

4

6

8

0

10

(N2)∞=1170[µM] (N2)∞=1.170 [µM]

0

2

4

6

t

8

10

t

C 0.016

(N2)∞=1170[µM] (N2)∞=1.170[µM]

0.014

0.012

u

0.01

0.008

0.006

0.004

0.002

0

2

4

6

8

10

t

Figure 12: Time course of N O • (N3 (t)) concentration in the microaerobic, (N1 )∞ = 15[µM ], regime for a thin biofilm with initial radius r = 0.04 [mm] under limited and unlimited nutrient concentration. A Robin boundary condition is imposed on DPTA in order to add it continuously to the system. Panel (A) shows the result for N O • (N3 (t)) and panel B displays the result for Hmp (H(t) produced and panel (C) shows the result for biofilm concentration u(t). The results are for (N4 )0 = 50[µM ]. 559

Biofilms are characterized by their complex heterogeneous structures and 37

t = 0.5[h]

t = 2[h]

t = 4[h]

t = 8[h]

Figure 13: Snapshot of Hmp(H) and oxygen concentration at t = 0.5, 2, 4, 8[h]. The color coding refers to the Hmp concentration and greyscal isolines show the oxygen concentration. Oxygen bulk concentration is (N1 )∞ = 15[µM ] and glucose bulk concentration is set at (N2 )∞ = 1.17[µM ].

38

560

substrate gradients. The outer surface has access to more nutrient while the

561

inner region may experience severe nutrient shortage, especially in an envi-

562

ronment with limited nutrient which results in the development of biofilms

563

with heterogeneous structures. To illustrate how the availability of nutrient

564

supply affects the spatial distribution of Hmp within the biofilm, two dimen-

565

sional visualisations are represented in Figure 13. The color coding refers to

566

the relative fraction of Hmp and greyscale isolines show the oxygen concen-

567

tration. We assume that the initial radius of the biofilm is r = 0.2 [mm]

568

and the oxygen bulk concentration is set at (N1 )∞ = 15[µM ]. The glucose

569

bulk concentration is chosen to be (N2 )∞ = 1.17[µM ]. Nutrient limitation

570

triggers the development of biofilms with heterogeneous distribution of Hmp,

571

c.f. Figure 13. The outer surface of the biofilm is mostly composed of Hmp

572

while in the inner region the density of Hmp is very low. Under nutrient lim-

573

ited conditions, more oxygen is available in the biofilm which results in the

574

production of Hmp in the region with high N O•. Since all substrates diffuse

575

into the biofilm from the aqueous phase, N O• concentration is higher near

576

the biofilm/liquid interface, c.f Figure 14. Therefore, more Hmp is produced

577

in this region which prevents N O• from reaching the inner layers of biofilms

578

and thus the density of biomass species Hmp decreases in this area.

579

5

Discussion

580

The mathematical model developed in this work is a novel approach to simu-

581

late N O• dynamics in growing biofilms. An important aspect of the proposed 39

t = 0.5

t=2

t=4

t=8

Figure 14: Snapshot of N O• concentration at t = 0.5, 2, 4, 8. The color coding refers to the N O• concentration. Oxygen bulk concentration is (N1 )∞ = 15[µM ] and glucose bulk concentration is set at (N2 )∞ = 1.17[µM ].

40

582

model is consideration of the complexity of bacterial biofilm spatial struc-

583

ture that improves the accuracy of N O• diffusion as well as other molecules

584

that impact biofilm development. The complexity of bacterial biofilm mi-

585

croanatomy has been recently illustrated by Serra et al. [44], who uncovered

586

spatial variation in the structure of E. coli biofilms and region-specific pat-

587

terns of cell differentiation. Due to the broad reactivity of N O• and system-

588

level perturbation that occurs in bacteria after N O• exposure, computation-

589

ally predicting the global response of a bacterial biofilm to N O• is highly

590

complex. To handle this difficulty, we restricted our model to incorporate a

591

set of the essential factors based on previous experimental studies. Of course,

592

there may be many further variations that should be taken into considera-

593

tion to modify the proposed model. Nevertheless, even in the current basic

594

model, we observe many interesting behaviors which estimate consequences

595

of N O• treatment. The conclusions reached by N O• modelling are in qual-

596

itative agreement with previous studies that have examined N O• dynamics

597

in planktonic bacteria [17, 41, 42], for instance, predicting the oscillatory be-

598

haviour of N O• under microaerobic conditions [41]. In addition, our model

599

may be useful to inform strategies for controlling unwanted biofilm forma-

600

tion in medical scenarios (e.g., bacterial infection, surgical wound healing),

601

as well as in industrial settings (biofouling in industrial and drinking water

602

systems).

41

603

N O• in a biological and medical context

604

Reactive nitrogen species form a portion of the human immune response to

605

persistent infections. Macrophage-generated nitric oxide acts to kill bacteria

606

through a variety of means. These include inhibition of vital processes and

607

damaging DNA [27]. In the current work, we specifically studied the detoxifi-

608

cation system and explored how to impair the enzymatic N O• detoxification

609

which consequently makes the biofilms more vulnerable against N O•. Ex-

610

isting work has also demonstrated the in vitro effectiveness of nitric oxide in

611

eradicating established biofilms, as well as rendering them more vulnerable

612

to antimicrobial agents [36]. In pathogenic biofilms formed by Pseudomonas

613

aeruginosa and other species, N O• is known to induce the dispersal of cells

614

in biofilms by lowering levels of the key second messenger molecule, bis-

615

(3 − 5)-cyclic dimeric GMP (c-di-GMP), which regulates the bacterial tran-

616

sition from a planktonic to sessile (biofilm) state [5]. Because of the effect of

617

N O• on biofilm dispersal, there has been active research on the use of N O•-

618

based antibiofilm strategies, including its use in combination with synthetic

619

antibiotics and native antimicrobials. Without nitric oxide exposure, most

620

antibiotics are completely ineffective at reducing established biofilms [36].

621

At present, the primary tools used to generate nitric oxide for medical

622

purposes are nitric oxide donor drugs. While effective in their production, the

623

ability of these substances to effectively target only the biofilm of interest is

624

limited: Nitric oxide disperses quickly, is short-lived, and highly reactive [30].

625

As the nitric oxide does not distinguish friend from foe, there is no real 42

626

option to provide a dose at a systemic level without risk of severe harm to

627

the patient [27]. At present, a major research trend is inclusion of these

628

substances in polymers in association with medical implants or at the site of

629

wounds, to allow for steady release. These substances are expensive, however,

630

and usually unstable, making their deployment highly selective [23].

631

Mathematical models of nitric oxide in biofilms provide an opportunity

632

to calculate the effective dose necessary to eliminate biofilms without ad-

633

ditional risk to the patient and with reduced cost for the treatment. This

634

ameliorates secondary harmful side effects of the treatment and its associated

635

risks. Accurate dose prediction could also minimize the amount of antibiotics

636

necessary in such treatments, lowering the likelihood of the development of

637

antibiotic resistance by pathogenic bacteria. Based on the obtained numeri-

638

cal results from our proposed model some strategies that impair the activity

639

of virulence factors and maximize the inhibitory effects of N O• against bac-

640

terial biofilms are increasing the initial concentration of electron donor or its

641

dissociation rate, treating the biofilm colony by N O• at the early stage of

642

its development, and implementing nutrient deprivation stress, see Figures 9

643

and 12.

644

Future directions and considerations

645

The results of this study show a link between biofilm resilience, oxygen avail-

646

ability, and nutrient availability. Already, this suggests a direction for ex-

647

perimental exploration here; testing whether patient starvation or increased 43

648

oxygenation might increase the effectiveness of nitric oxide donor drugs. Fu-

649

ture computational explorations will rely on further developments of the

650

model. It is known from prior experimental investigations that biofilms are

651

not uniform in oxygenation, with deeper layers tending to be more oxygen-

652

starved [1]. This generates a gradient of oxygen, producing a range of oxy-

653

gen environments, ranging from approximating the background to completely

654

anoxic states. This suggests that, when any amount of oxygen is present, a

655

biofilm is vulnerable to attack: Some microaerobic "zone" will exist between

656

the media interface and the anoxic zone within the biofilm. This opens a new

657

avenue of investigation into whether biofilms may actively generate anaero-

658

bic environments to diminish the effectiveness of oxidative attack. Consistent

659

with this, studies have shown that P. aeruginosa PA01 produces a highly ex-

660

clusive mucus layer that may act as a defense against attack by macrophages.

661

This mucus layer is absent in anaerobic conditions [43].

662

There are also species-specific modifications to this model that inevitably

663

arise. However, these modifications are likely only to be worth the inves-

664

tigative effort for specific species. These include, for example, variations in

665

composition, density, and flow-rate, and the effects of biofilm age on each, all

666

of which would add additional variables to the effectiveness and penetration

667

of the nitric oxide. Work by de Beer and colleagues showed that the age

668

of biofilms affects their density and oxygen gradient, but also demonstrated

669

that species-specific variations to the overall structure - pores, etc. - greatly

670

alter this as well, no doubt affecting the rate of dispersal of N O2 through the 44

671

biofilm in the process [1].

672

There are also considerations of nitric oxide interactions as a part of this

673

process. In studies investigating the action of macrophages, it will be insuf-

674

ficient only to consider the effects of nitric oxide. Nitric oxide is generated

675

along with other radical species, including superoxides, which can react with

676

nitric oxide to form peroxynitrite [3]. This substance is exceedingly harm-

677

ful to bacteria, as it can induce strand breaks in their DNA [3]. However,

678

many bacteria also release eDNA as a part of their biofilms - some even have

679

DNA as a primary constituent of the ECM [11,46]. This is particularly well-

680

known in S. aureus, a major human pathogen [46]. It is plausible that this

681

eDNA may act as another kind of barrier defense; depleting the nitric oxide,

682

superoxide, and peroxynitrite before it has a chance to act on the bacteria

683

within the biofilm. Future mathematical models that include these additional

684

factors may have increased accuracy and potential in predicting macrophage

685

behaviour. Ultimately, however, the model presented here provides a starting

686

framework, on which more detailed and parameter-rich models of molecule

687

dynamics within biofilms can be built.

688

6

Conclusion

689

N O• is a potent antimicrobial which is able to diffuse through cellular mem-

690

branes and activate diverse effects. Its role in defending against infection is

691

critical. At high concentrations, N O• binds DNA, proteins and lipids and

692

as a results kills or inhibits target pathogens. The broad reactivity of N O• 45

693

makes it difficult to predict its biological effect on bacteria under different

694

environmental conditions. Our primary objective in this paper has been to

695

introduce a mathematical model that can simulate the dynamics of N O• in

696

a growing biofilm. To validate the model we studied the N O• stress response

697

mechanism under aerobic, microaerobic and anaerobic conditions. The ob-

698

tained results are in good agreement qualitatively with experimental results

699

for well mixed cultures in the literature. We investigated the effects of initial

700

concentration levels of N O• donor and its half-life time, influence of glucose

701

bulk concentration levels as the only carbon source, and initial size of biofilm

702

colony on inhibiting the virulence factors. The main findings are:

703

• Oscillatory behavior of N O• concentration under a microaerobic regime

704

which increases the clearance time of N O• is a temporary phenomenon

705

and does not generate gradients in N O• concentration levels within the

706

biofilm.

707

708

• In order to impair the system of N O• enzymatic detoxification, large amounts of N O• donor with short half-lives should be used.

709

• The initial size of a biofilm colony affects the N O• detoxification. By

710

decreasing the initial size of biofilm colonies, less N O• can diffuse into

711

the biofilm and its consumption by reducing agents also decreases.

712

Compared with the reduction in N O• detoxification by Hmp as an

713

aerobic reducing agent, the reduction in the flux of N O• is less pro-

714

nounced. This indicates that targeting biofilms at the early stage of 46

715

their growth can increase N O• level which eradicates the biofilms more

716

efficiently. The oscillatory behavior of N O• which manifests under mi-

717

croaerobic conditions disappears in cultures with small initial size due

718

to the formation of small amounts of Hmp that can not compete with

719

other oxygen dependent activities which consume oxygen.

720

• Nutrient deprivation stress can maximize the inhibitory effect of N O•

721

against bacterial biofilms. The effect of nutrient shortage on removing

722

biofilms by N O• is more pronounced under microaerobic conditions and

723

if N O• donor is added continuously (using Robin boundary condition)

724

at the early stage of biofilm growth. Under these conditions, production

725

of Hmp is decreased which maximizes the biofilm removal by N O•.

726

The obtained results highlight the role of oxygen, biofilm size and the

727

amount of N O• in the system, on the biofilm eradication by N O•.

728

• For thick biofilms under microaerobic regimes, nutrient limitation pre-

729

vents N O• from diffusing into the inner region of the biofilm and leads

730

to the layering phenomenon such that the outer layers of the biofilm

731

consist of Hmp while in the inner regions the density of Hmp is neg-

732

ligible. This emphasizes the importance of spatio-temporal modeling

733

to study the dynamics of N O• and its detoxification system. It also

734

reveals the reason of resistance of thick biofilms against removal. The

735

inner layers do not have access to enough N O• and Hmp detoxifies

736

N O• at the outer surface of biofilms which consequently impair biofilm 47

737

eradication.

738

Acknowledgement. This research was supported in parts by the Natu-

739

ral Science and Engineering Research Council of Canada (NSERC) with a

740

Discovery Grant awarded to SS. We also thank the Fields institute and Uni-

741

versity of Waterloo for providing financial support for MG as a postdoctoral

742

fellow.

743

References

744

[1] Beer, D. D., Stoodley, P., Roe, F., Lewandowski, Z. (1994).

745

Effects of biofilm structures on oxygen distribution and mass

746

transport. Biotechnology and Bioengineering. 43(11),

747

doi:10.1002/bit.260431118.

1131-1138.

748

[2] Bowman, L.A., McLean, S., Poole, R.K., Fukuto, J.M. (2011.) The Di-

749

versity of Microbial Responses to Nitric Oxide and Agents of Nitrosative

750

Stress: Close Cousins but Not Identical Twins. Adv. Microb. Physiol. 59,

751

135-219.

752

[3] Burney, S., Caulfield, J. L., Niles, J. C., Wishnok, J. S., Tannenbaum,

753

S. R. (1999). The chemistry of DNA damage from nitric oxide and per-

754

oxynitrite. Mutation Research/Fundamental and Molecular Mechanisms

755

of Mutagenesis. 424(1-2), 37-49. doi:10.1016/s0027-5107(99)00006-8. 48

756

[4] Cegelski, L., Marshall,G.R., Eldridge, G.R., Hultgren, S.J. (2008). The

757

biology and future prospects of antivirulence therapies. Nature Review

758

Microbiology. 6(1), 17-27.

759

[5] Cutruzzolá, F., Frankenberg-Dinkel, N. (2015) Origin and Impact of Ni-

760

tric Oxide in Pseudomonas aeruginosa Biofilms. Journal of Bacteriology.

761

198(1): 55-65. DOI: 10.1128/JB.00371-15.

762

[6] Eberl, H. J., Parker, D. F., Van Loosdrecht, C. M. (2001). A new deter-

763

ministic spatio-temporal continuum model for biofilm development. J.

764

of Theoretical Medicine. 3, 161-175.

765

[7] Eberl, H.J., Demaret, L. (2007). A finite difference scheme for a degen-

766

erated diffusion equation arising in microbial ecology. El. J. Diff. Equs.

767

CS 15, 77-95.

768

[8] Eberl, H. J., Collinson, S. (2009). A modelling and simulation study of

769

siderophore mediated antagonsim in dual-species biofilms. Theor. Biol.

770

Med. Mod. 6:30.

771

[9] Efendiev, M.A., Zelik, S.V., Eberl, H.J. (2009). Existence and longtime

772

behaviour of a biofilm model. Comm. Pure and Appl. Analysis. 8(2),

773

509-531.

774

775

[10] Fernebro, J. (2011). Fighting bacterial infections-future treatment options. Drug Resist Updates. 14(2), 125-139. 49

776

777

[11] Flemming, H.-C., Wingender, J. (2010). The biofilm matrix. Nature Reviews Microbiology. 8(9), 623-633. doi: 10.1038/nrmicro2415.

778

[12] Frederick, M., Kuttler, C., Hense, B.A., Müller, J., Eberl, H.J. (2010). A

779

mathematical model of quorum sensing in patchy biofilm communities

780

with slow background flow. Can. Appl. Math. Quart. 18(3), 267-298.

781

[13] Ghasemi, M., Eberl, H. J. (2017). Time adaptive numerical solution of

782

a highly degenerate diffusion-reaction biofilm model based on regulari-

783

sation. J. Sci. Comp. 74, 1060-1090.

784

[14] Ghasemi, M., Eberl, H. J. (2017). Extension of a Regularization Based

785

Time-adaptive Numerical Method for a Degenerate Diffusion-Reaction

786

Biofilm Growth Model to Systems Involving Quorum Sensing. Proc.

787

Comp. Sci. 108, 1893-1902.

788

[15] Ghasemi, M., Hense, B.A., Eberl, H. J, Kuttler, C. (2018). Simulation-

789

Based Exploration of Quorum Sensing Triggered Resistance of Biofilms

790

to Antibiotics. Bull. Math. Biol. 80, 1736-1775.

791

[16] Ghasemi, M., Sonner, S., Eberl, H. J. (2018). Time adaptive numerical

792

solution of a highly nonlinear degenerate cross-diffusion system arising

793

in multi-species biofilm modeling Eur. J. Appl. Math. 29, 1035-1061.

794

[17] Gowers, G.F., Robinson, J.L., Brynildsen. M.P. (2016). Starved Es-

795

cherichia coli preserve reducing power under nitric oxide stress. Bio-

796

chemical and Biophysical Research Communications. 476, 29-32. 50

797

[18] Hall-Stoodley, L., Costerton, J.W., Stoodley, P. (2004). Bacterial

798

biofilms: from the natural environment to infectious diseases. Nature

799

Reviews Microbiology. 2(2), 95-108.

800

[19] Hu, T.M., Hayton, W.L., Mallery, S.R. (2006). Kinetic Modeling of

801

Nitric-Oxide- Associated Reaction Network.Pharm. Res. 23, 1702-1711.

802

[20] Khassehkhan, H., Hillen, T., Eberl, H.J. (2009) A non-linear master

803

equation for a degenerate diffusion model of biofilm growth. LNCS. 5544,

804

735-744.

805

[21] Lancaster, J.R, Jr. (2006) Nitroxidative, Nitrosative, and Nitrative

806

Stress: Kinetic Predictions of Reactive Nitrogen Species Chemistry Un-

807

der Biological Conditions. Chem. Res. Toxicol. 19, 1160-1174.

808

809

[22] Lear, G., Lewis, G.D. (2012). Microbial Biofilms: Current Research and Applications. Caister Academic Press. ISBN 978-1-904455-96-7.

810

[23] Lee, W. H., Ren, H., Wu, J., Novak, O., Brown, R. B., Xi, C., Meyerhoff,

811

M. E. (2016). Electrochemically Modulated Nitric Oxide Release From

812

Flexible Silicone Rubber Patch: Antimicrobial Activity For Potential

813

Wound Healing Applications. ACS Biomaterials Science and Engineer-

814

ing. 2(9), 1432-1435. doi:10.1021/acsbiomaterials.6b00360.

815

[24] Lewis, R.S., Deen, W.M. (1994). Kinetics of the reaction of Nitric Oxide

816

with oxygen in aqueous solutions. Chem. Res. Toxicol. 7, 568-574. 51

817

[25] Lewis, R.S., Tamir, S., Tannenbaum, S.R., Deen, W.M. (1995). Kinetic

818

Analysis of the Fate of Nitric Oxide Synthesized by Macrophages in

819

Vitro.J. Biol. Chem. 270, 29350-29355.

820

[26] Lim, C.H., Dedon, P.C., Deen, W.A. (2008). Kinetic Analysis of Intracel-

821

lular Concentrations of Reactive Nitrogen Species. Chem. Res. Toxicol.

822

21, 2134-2147.

823

[27] MacMicking, J., Xie, Q., Nathan, C. (1997). Nitric Oxide And

824

Macrophage Function. Annual Review of Immunology. 15(1), 323-350.

825

doi:10.1146/annurev.immunol.15.1.323.

826

[28] Mason, M.G., Holladay, R.S., Nicholls, P., Shepherd, M., Cooper, C.E.

827

(2008). A quantitative approach to Nitric Oxide inhibition of terminal

828

oxidases of the respiratory chain. Methods Enzymol. 437, 135-159.

829

[29] Mason, M.G., Shepherd, M., Nicholls, P.,Dobbin, P.S. (2009). Cy-

830

tochrome bd confers nitric oxide resistance to Escherichia coli. Nat.

831

Chem. Biol. 5, 94-96.

832

[30] Miller, M. R., Megson, I. L. (2009). Recent developments in nitric

833

oxide donor drugs. British Journal of Pharmacology, 151(3), 305-321.

834

doi:10.1038/sj.bjp.0707224.

835

[31] Piotrowska, M.J., BartLomiejczyk, A., Bodnar, M. (2018). Mathemati-

836

cal analysis of a generalised p53-Mdm2 protein gene expression model.

837

Applied Mathematics and Computation. 328, 26-44. 52

838

839

[32] Ngamsaad, W., Sunatai, S. (2016). Mechanically-driven spreading of bacterial populations. Comm. Nonlin. Sci. Num. Sim. 35, 88-96.

840

[33] Poole, R.K., Cook, G.M. (2000). Redundancy of aerobic respiratory

841

chains in bacteria? routes, reasons and regulation. Adv. Microb. Physiol.

842

43, 165-224.

843

844

845

846

[34] Rang, J. Improved traditional Rosenbrock-Wanner methods for stiff ODEs and DAEs. Institute of scientific computing, Germany, 2013. [35] Reiter, T.A. (2006). N O∗ chemistry: a diversity of targets in the cell. Redox. Rep. 11, 194-206.

847

[36] Ren, H., Wu, J., Colletta, A., Meyerhoff, M. E., Xi, C. (2016). Efficient

848

Eradication of Mature Pseudomonas aeruginosa Biofilm via Controlled

849

Delivery of Nitric Oxide Combined with Antimicrobial Peptide and An-

850

tibiotics. Frontiers in Microbiology. 7. doi:10.3389/fmicb.2016.01260

851

[37] Robinson, J.L., Brynildsen, M.P. (2013) A Kinetic Platform to Deter-

852

mine the Fate of Nitric Oxide in Escherichia coli. PLoS Comput Biol

853

9(5): e1003049. doi:10.1371/journal.pcbi.1003049.

854

[38] Robinson, J.L., Miller, R.V., Brynildsen, M.P. (2014) Model-driven

855

identification of dosing regimens that maximize the antimicrobial ac-

856

tivity of nitric oxide. Metabolic Engineering Communication 1, 12-18. 53

857

[39] Robinson, J.L., Adolfsen, K.J, Brynildsen, M.P. (2014) Diciphering ni-

858

tric oxide stress in bacteria with quantitative modeling Current Opinion

859

in Microbiology 19, 16-24.

860

[40] Robinson, J.L., Brynildsen, M.P. Ensemble Modeling Enables Quanti-

861

tative Exploration of Bacterial Nitric Oxide Stress Networks. Stress and

862

Environmental Control of Gene Expression and Adaptation in Bacteria.

863

Chapter 17.6, 1009-1014.

864

[41] Robinson, J.L., Brynildsen, M.P. (2016). Discovery and dissection of

865

metabolic oscillations in the microaerobic nitric oxide response network

866

of Escherichia coli. PNAS 1757-1766

867

[42] Robinson, J.L., Jaslove, J.M., Murawski, A.M., Fazen, C.H., Brynildsen,

868

M.P. (2017). An integrated network analysis reveals that nitric oxide

869

reductase prevents metabolic cycling of nitric oxide by Pseudomonas

870

aeruginosa. Metabolic Engineering 41, 67-81.

871

[43] Sabra, W., Kim, E.J., Zeng, A.P. (2002). Physiological responses of

872

Pseudomonas aeruginosa PAO1 to oxidative stress in controlled mi-

873

croaerobic and aerobic cultures. Microbiology. 148, 3195-202.

874

[44] Serra, D.O., Richter, A.M., Klauck, G., Mika, F., Hengge, R. (2013).

875

Microanatomy at Cellular Resolution and Spatial Order of Physiolog-

876

ical Differentiation in a Bacterial Biofilm. mBio. 4(2), e00103-13. doi:

877

10.1128/mBio.00103-13 54

878

[45] Shepherd, M., Achard, m.e.s, Idris, A.,Totsika, M., Phan, M., Pe-

879

ters, K.M., Sarkar, S., Ribeiro1, C.A., Holyoake1, L.V., Ladakis1, D.,

880

Ulett, G.C., Sweet, M.J., Poole, R.K., McEwan, A.G., Schembri, M.A.

881

(2016). The cytochrome bd-I respiratory oxidase augments survival of

882

multidrug-resistant Escherichia coli during infection. Scientific Reports.

883

6:35285, DOI: 10.1038/srep35285.

884

[46] Sugimoto, S., Sato, F., Miyakawa, R., Chiba, A., Onodera, S., Hori, S.,

885

Mizunoe, Y. (2018). Broad impact of extracellular DNA on biofilm for-

886

mation by clinically isolated Methicillin-resistant and -sensitive strains

887

of Staphylococcus aureus. Scientific Reports. 8(1). doi: 10.1038/s41598-

888

018-20485-z.

889

[47] Toledo, J.C. Jr., Augusto, O. (2012). Connecting the chemical and bio-

890

logical properties of nitric oxide. Chem. Res. Toxicol. 25, 975-989.

891

[48] Wanner, O., Eberl, H. J., Van Loosdrecht, M. C. M., Morgenroth,

892

E., Noguera, D. R., Picioreanu, C, Rittmann, BE. (2006).Mathemati-

893

cal Modelling of Biofilms, IWA Publishing.

894

895

[49] Watnick, P., Kolter, R. (2000). Biofilm – City of Microbes (Minireview). J. Bacteriol. 182(10), 2675-2679.

55