kinetic_formalism.Rnw

   1 \documentclass[english,12pt]{article}
   2 \usepackage{fancyhdr}
   3 %\usepackage[pdftex]{graphicx}
   4 \usepackage{graphicx}
   5 \usepackage[bf]{caption2}
   6 \usepackage{rotating}
   7 \usepackage{multirow}
   8 \usepackage{textcomp}
   9 \usepackage{mathrsfs}
  10 \usepackage{amssymb}
  11 \usepackage{setspace}
  12 \usepackage{txfonts}
  13 \usepackage[light,all]{draftcopy}
  14 \usepackage{fancyref}
  15 \usepackage[hyperfigures,backref,bookmarks,colorlinks]{hyperref}
  16 \usepackage[sectionbib,sort&compress,square,numbers]{natbib}
  17 \usepackage[margin,inline,draft]{fixme}
  18 \usepackage[x11names,svgnames]{xcolor}
  19 \usepackage{texshade}
  20 \newenvironment{narrow}[2]{%
  21   \begin{list}{}{%
  22       \setlength{\topsep}{0pt}%
  23       \setlength{\leftmargin}{#1}%
  24       \setlength{\rightmargin}{#2}%
  25       \setlength{\listparindent}{\parindent}%
  26       \setlength{\itemindent}{\parindent}%
  27       \setlength{\parsep}{\parskip}}%
  28   \item[]}{\end{list}}
  29 \newenvironment{paperquote}{%
  30   \begin{quote}%
  31      \it
  32   }%
  33   {\end{quote}}
  34 \renewcommand{\textfraction}{0.15}
  35 \renewcommand{\topfraction}{0.85}
  36 \renewcommand{\bottomfraction}{0.65}
  37 \renewcommand{\floatpagefraction}{0.60}
  38 %\renewcommand{\baselinestretch}{1.8}
  39 \newenvironment{enumerate*}%
  40   {\begin{enumerate}%
  41     \setlength{\itemsep}{0pt}%
  42     \setlength{\parskip}{0pt}}%
  43   {\end{enumerate}}
  44 \newenvironment{itemize*}%
  45   {\begin{itemize}%
  46     \setlength{\itemsep}{0pt}%
  47     \setlength{\parskip}{0pt}}%
  48   {\end{itemize}}
  49 \oddsidemargin 0.0in
  50 \textwidth 6.5in
  51 \raggedbottom
  52 \clubpenalty = 10000
  53 \widowpenalty = 10000
  54 \pagestyle{fancy}
  55 \author{Don Armstrong}
  56 \title{OOL Kinetic Formalisms}
  57 %\date{}
  58 \onehalfspacing
  59 \begin{document}
  60 %\maketitle
  61
  62 <<results=hide,echo=FALSE>>=
  63 require(lattice)
  64 require(grid)
  65 # R in cal / mol K
  66 to.kcal <- function(k,temp=300) {
  67   gasconst <- 1.985
  68   return(-gasconst*temp*log(k)/1000)
  69 }
  70 @
  71
  72 \section{State Equation}
  73 % double check this with the bits in the paper
  74
  75 Given a base forward kinetic parameter for the $i$th specie $k_{fi}$
  76 (which is dependent on lipid type, that is PC, PE, PS, etc.), an
  77 adjustment parameter $k_{fi\mathrm{adj}}$ based on the vesicle and the
  78 specific specie (length, unsaturation, etc.) (see~\fref{eq:kf_adj}),
  79 the molar concentration of monomer of the $i$th specie
  80 $\left[C_{i_\mathrm{monomer}}\right]$, the surface area of the vesicle
  81 $S_\mathrm{ves}$, the base backwards kinetic parameter for the $i$th
  82 specie $k_{bi}$ which is also dependent on lipid type, its adjustment
  83 parameter $k_{bi\mathrm{adj}}$ (see~\fref{eq:kb_adj}), and the molar
  84 concentration of the $i$th specie in the vesicle $C_{i_\mathrm{ves}}$,
  85 the change in concentration of the $i$th specie in the vesicle per
  86 change in time $\frac{d C_{i_\mathrm{ves}}}{dt}$ can be calculated:
  87
  88 \begin{equation}
  89   \frac{d C_{i_\mathrm{ves}}}{dt} = k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]S_\mathrm{ves} -
  90   k_{bi}k_{bi\mathrm{adj}}C_{i_\mathrm{ves}}
  91   \label{eq:state}
  92 \end{equation}
  93
  94 For $k_{fi}k_{fi\mathrm{adj}}\left[C_{i_\mathrm{monomer}}\right]$,
  95 $k_{fi}$ has units of $\frac{\mathrm{m}}{\mathrm{s}}$,
  96 $k_{fi\mathrm{adj}}$ and $k_{bi\mathrm{adj}}$ are unitless,
  97 concentration is in units of $\frac{\mathrm{n}}{\mathrm{L}}$, surface
  98 area is in units of $\mathrm{m}^2$, $k_{bi}$ has units of
  99 $\frac{1}{\mathrm{s}}$ and $C_{i_\mathrm{ves}}$ has units of
 100 $\mathrm{n}$, Thus, we have
 101
 102 \begin{equation}
 103   \frac{\mathrm{n}}{\mathrm{s}} = \frac{\mathrm{m}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \mathrm{m}^2 \frac{1000\mathrm{L}}{\mathrm{m}^3} -
 104   \frac{1}{\mathrm{s}} \mathrm{n}
 105   =
 106   \frac{\mathrm{m^3}}{\mathrm{s}} \frac{\mathrm{n}}{\mathrm{L}} \frac{1000\mathrm{L}}{\mathrm{m}^3} - \frac{\mathrm{n}}{\mathrm{s}}=
 107   \frac{\mathrm{n}}{\mathrm{s}} = 1000 \frac{\mathrm{n}}{\mathrm{s}} - \frac{\mathrm{n}}{\mathrm{s}}
 108   \label{eq:state_units}
 109 \end{equation}
 110
 111 The 1000 isn't in \fref{eq:state} above, because it is unit-dependent.
 112
 113 \subsection{Forward adjustments ($k_{fi\mathrm{adj}}$)}
 114
 115 The forward rate constant adjustment, $k_{fi\mathrm{adj}}$ takes into
 116 account unsaturation ($un_f$), charge ($ch_f$), curvature ($cu_f$),
 117 length ($l_f$), and complex formation ($CF1_f$), each of which are
 118 modified depending on the specific specie and the vesicle into which
 119 the specie is entering.
 120
 121 \begin{equation}
 122   k_{fi\mathrm{adj}} = un_f \cdot ch_f \cdot cu_f \cdot l_f \cdot CF1_f
 123   \label{eq:kf_adj}
 124 \end{equation}
 125
 126 \newpage
 127 \subsubsection{Unsaturation Forward}
 128
 129 In order for a lipid to be inserted into a membrane, a void has to be
 130 formed for it to fill. Voids can be generated by the combination of
 131 unsaturated and saturated lipids forming herterogeneous domains. Void
 132 formation is increased when the unsaturation of lipids in the vesicle
 133 is widely distributed; in other words, the insertion of lipids into
 134 the membrane is greater when the standard deviation of the
 135 unsaturation is larger. Assuming that an increase in width of the
 136 distribution linearly decreases the free energy of activation, the
 137 $un_f$ parameter must follow
 138 $x^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}$ where $x > 1$, so a
 139 convenient starting base for $x$ is 2:
 140
 141 \begin{equation}
 142   un_f = 2^{\mathrm{stdev}\left(un_\mathrm{ves}\right)}
 143   \label{eq:unsaturation_forward}
 144 \end{equation}
 145
 146 \setkeys{Gin}{width=3.2in}
 147 <<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
 148 curve(2^x,from=0,to=sd(c(0,4)),
 149       main="Unsaturation Forward",
 150       xlab="Standard Deviation of Unsaturation of Vesicle",
 151       ylab="Unsaturation Forward Adjustment")
 152 @
 153 <<fig=TRUE,echo=FALSE,results=hide,width=5,height=5>>=
 154 curve(to.kcal(2^x),from=0,to=sd(c(0,4)),
 155       main="Unsaturation forward",
 156       xlab="Standard Deviation of Unsaturation of Vesicle",
 157       ylab="Unsaturation Forward (kcal/mol)")
 158 @
 159
 160
 161 \newpage
 162 \subsubsection{Charge Forward}
 163
 164 A charged lipid such as PS approaching a vesicle with an average
 165 charge of the same sign will experience repulsion, whereas those with
 166 different signs will experience attraction, the degree of which is
 167 dependent upon the charge of the monomer and the average charge of the
 168 vesicle. If either the vesicle or the monomer has no charge, there
 169 should be no effect of charge upon the rate. This leads us to the
 170 following equation, $x^{-\left<ch_v\right> ch_m}$, where
 171 $\left<ch_v\right>$ is the average charge of the vesicle, and $ch_m$
 172 is the charge of the monomer. If either $\left<ch_v\right>$ or $ch_m$
 173 is 0, the adjustment parameter is 1 (no change), whereas it decreases
 174 if both are positive or negative, as the product of two real numbers
 175 with the same sign is always positive. A convenient base for $x$ is
 176 60, resulting in the following equation:
 177
 178
 179 \begin{equation}
 180   ch_f = 60^{-\left<{ch}_v\right> {ch}_m}
 181   \label{eq:charge_forward}
 182 \end{equation}
 183
 184 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 185 x <- seq(-1,0,length.out=20)
 186 y <- seq(-1,0,length.out=20)
 187 grid <- expand.grid(x=x,y=y)
 188 grid$z <- as.vector(60^(-outer(x,y)))
 189 print(wireframe(z~x*y,grid,cuts=50,
 190                 drape=TRUE,
 191                 scales=list(arrows=FALSE),
 192                 main="Charge Forward",
 193                 xlab=list("Average Vesicle Charge",rot=30),
 194                 ylab=list("Component Charge",rot=-35),
 195                 zlab=list("Charge Forward",rot=93)))
 196 rm(x,y,grid)
 197 @
 198 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 199 x <- seq(-1,0,length.out=20)
 200 y <- seq(-1,0,length.out=20)
 201 grid <- expand.grid(x=x,y=y)
 202 grid$z <- as.vector(to.kcal(60^(-outer(x,y))))
 203 print(wireframe(z~x*y,grid,cuts=50,
 204                 drape=TRUE,
 205                 scales=list(arrows=FALSE),
 206                 main="Charge Forward (kcal/mol)",
 207                 xlab=list("Average Vesicle Charge",rot=30),
 208                 ylab=list("Component Charge",rot=-35),
 209                 zlab=list("Charge Forward (kcal/mol)",rot=93)))
 210 rm(x,y,grid)
 211 @
 212
 213
 214 \newpage
 215 \subsubsection{Curvature Forward}
 216
 217 Curvature is a measure of the intrinsic propensity of specific lipids
 218 to form micelles (positive curvature), inverted micelles (negative
 219 curvature), or planar sheets (zero curvature). In this formalism,
 220 curvature is measured as the ratio of the size of the head to that of
 221 the base, so negative curvature is bounded by $(0,1)$, zero curvature
 222 is 1, and positive curvature is bounded by $(1,\infty)$. The curvature
 223 can be transformed into the typical postive/negative mapping using
 224 $\log$, which has the additional property of making the range of
 225 positive and negative curvature equal, and distributed about 0.
 226
 227 As in the case of unsaturation, void formation is increased by the
 228 presence of lipids with mismatched curvature. Thus, a larger
 229 distribution of curvature in the vesicle increases the rate of lipid
 230 insertion into the vesicle. However, a species with curvature $e^-1$
 231 will cancel out a species with curvature $e$, so we have to log
 232 transform (turning these into -1 and 1), then take the absolute value
 233 (1 and 1), and finally measure the width of the distribution. Thus, by
 234 using the log transform to make the range of the lipid curvature equal
 235 between positive and negative, taking the absolute value to cancel out
 236 exactly mismatched curvatures, we come to an equation with the shape
 237 $x^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}$, and a
 238 convenient base for $x$ is 10, yielding:
 239
 240 {\color{red} Shouldn't a vesicle of -1,1,0 have the same activation
 241   energy as a vesicle of 0,0,0? It doesn't currently.}
 242
 243
 244 \begin{equation}
 245   cu_f = 10^{\mathrm{stdev}\left|\log cu_\mathrm{vesicle}\right|}
 246   \label{eq:curvature_forward}
 247 \end{equation}
 248
 249
 250 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
 251 curve(10^x,from=0,to=max(c(sd(abs(log(c(0.8,1.33)))),
 252                     sd(abs(log(c(1,1.33)))),
 253                     sd(abs(log(c(0.8,1)))))),
 254       main="Curvature forward",
 255       xlab="Standard Deviation of Absolute value of the Log of the Curvature of Vesicle",
 256       ylab="Curvature Forward Adjustment")
 257 @
 258 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
 259 curve(to.kcal(10^x),from=0,to=max(c(sd(abs(log(c(0.8,1.33)))),
 260                     sd(abs(log(c(1,1.33)))),
 261                     sd(abs(log(c(0.8,1)))))),
 262       main="Curvature forward",
 263       xlab="Standard Deviation of Absolute value of the Log of the Curvature of Vesicle",
 264       ylab="Curvature Forward Adjustment (kcal/mol)")
 265 @
 266
 267
 268 \newpage
 269 \subsubsection{Length Forward}
 270
 271 As in the case of unsaturation, void formation is easier when vesicles
 272 are made up of components of widely different lengths. Thus, when the
 273 width of the distribution of lengths is larger, the forward rate
 274 should be greater as well, leading us to an equation of the form
 275 $x^{\mathrm{stdev} l_\mathrm{ves}}$, where $\mathrm{stdev}
 276 l_\mathrm{ves}$ is the standard deviation of the length of the
 277 components of the vesicle, which has a maximum possible value of 8 and
 278 a minimum of 0 in this set of experiments. A convenient base for $x$
 279 is 2, leading to:
 280
 281 \begin{equation}
 282   l_f = 2^{\mathrm{stdev} l_\mathrm{ves}}
 283   \label{eq:length_forward}
 284 \end{equation}
 285
 286 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
 287 curve(2^x,from=0,to=sd(c(12,24)),
 288       main="Length forward",
 289       xlab="Standard Deviation of Length of Vesicle",
 290       ylab="Length Forward Adjustment")
 291 @
 292 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=5>>=
 293 curve(to.kcal(2^x),from=0,to=sd(c(12,24)),
 294       main="Length forward",
 295       xlab="Standard Deviation of Length of Vesicle",
 296       ylab="Length Forward Adjustment (kcal/mol)")
 297 @
 298
 299
 300 \subsubsection{Complex Formation}
 301 There is no contribution of complex formation to the forward reaction
 302 rate in the current formalism.
 303
 304 \begin{equation}
 305   CF1_f=1
 306   \label{eq:complex_formation_forward}
 307 \end{equation}
 308
 309 \subsection{Backward adjustments ($k_{bi\mathrm{adj}}$)}
 310
 311 Just as the forward rate constant adjustment $k_{fi\mathrm{adj}}$
 312 does, the backwards rate constant adjustment $k_{bi\mathrm{adj}}$
 313 takes into account unsaturation ($un_b$), charge ($ch_b$), curvature
 314 ($cu_b$), length ($l_b$), and complex formation ($CF1_b$), each of
 315 which are modified depending on the specific specie and the vesicle
 316 into which the specie is entering:
 317
 318
 319 \begin{equation}
 320   k_{bi\mathrm{adj}} = un_b \cdot ch_b \cdot cu_b \cdot l_b \cdot CF1_b
 321   \label{eq:kb_adj}
 322 \end{equation}
 323
 324 \subsubsection{Unsaturation Backward}
 325
 326 The reverse rate of the reaction
 327
 328 \begin{equation}
 329   un_b = 10^{\left|3.5^{-\left<un_\mathrm{ves}\right>}-3.5^{-un_\mathrm{monomer}}\right|}
 330   \label{eq:unsaturation_backward}
 331 \end{equation}
 332
 333 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 334 grid <- expand.grid(x=seq(0,4,length.out=20),
 335                     y=seq(0,4,length.out=20))
 336 grid$z <- 10^(abs(3.5^-grid$x-3.5^-grid$y))
 337 print(wireframe(z~x*y,grid,cuts=50,
 338           drape=TRUE,
 339           scales=list(arrows=FALSE),
 340           xlab=list("Average Vesicle Unsaturation",rot=30),
 341           ylab=list("Monomer Unsaturation",rot=-35),
 342           zlab=list("Unsaturation Backward",rot=93)))
 343 rm(grid)
 344 @
 345 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 346 grid <- expand.grid(x=seq(0,4,length.out=20),
 347                     y=seq(0,4,length.out=20))
 348 grid$z <- to.kcal(10^(abs(3.5^-grid$x-3.5^-grid$y)))
 349 print(wireframe(z~x*y,grid,cuts=50,
 350           drape=TRUE,
 351           scales=list(arrows=FALSE),
 352           xlab=list("Average Vesicle Unsaturation",rot=30),
 353           ylab=list("Monomer Unsaturation",rot=-35),
 354           zlab=list("Unsaturation Backward (kcal/mol)",rot=93)))
 355 rm(grid)
 356 @
 357
 358
 359 \newpage
 360 \subsubsection{Charge Backwards}
 361 \begin{equation}
 362   ch_b = 20^{\left<{ch}_v\right> {ch}_m}
 363   \label{eq:charge_backwards}
 364 \end{equation}
 365
 366 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 367 x <- seq(-1,0,length.out=20)
 368 y <- seq(-1,0,length.out=20)
 369 grid <- expand.grid(x=x,y=y)
 370 grid$z <- as.vector(20^(outer(x,y)))
 371 print(wireframe(z~x*y,grid,cuts=50,
 372           drape=TRUE,
 373           scales=list(arrows=FALSE),
 374           xlab=list("Average Vesicle Charge",rot=30),
 375           ylab=list("Component Charge",rot=-35),
 376           zlab=list("Charge Backwards",rot=93)))
 377 rm(x,y,grid)
 378 @
 379 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 380 x <- seq(-1,0,length.out=20)
 381 y <- seq(-1,0,length.out=20)
 382 grid <- expand.grid(x=x,y=y)
 383 grid$z <- to.kcal(as.vector(20^(outer(x,y))))
 384 print(wireframe(z~x*y,grid,cuts=50,
 385           drape=TRUE,
 386           scales=list(arrows=FALSE),
 387           xlab=list("Average Vesicle Charge",rot=30),
 388           ylab=list("Component Charge",rot=-35),
 389           zlab=list("Charge Backwards (kcal/mol)",rot=93)))
 390 rm(x,y,grid)
 391 @
 392
 393 \newpage
 394 \subsubsection{Curvature Backwards}
 395 \begin{equation}
 396   cu_f = 7^{1-\left(20\left(\log_{e} cu_\mathrm{vesicle}-\log_{e} cu_\mathrm{monomer}\right)^2+1\right)^{-1}}
 397   \label{eq:curvature_backwards}
 398 \end{equation}
 399
 400 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 401 grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
 402                     y=seq(0.8,1.33,length.out=20))
 403 grid$z <- 7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1))
 404 print(wireframe(z~x*y,grid,cuts=50,
 405           drape=TRUE,
 406           scales=list(arrows=FALSE),
 407           xlab=list("Vesicle Curvature",rot=30),
 408           ylab=list("Monomer Curvature",rot=-35),
 409           zlab=list("Curvature Backward",rot=93)))
 410 rm(grid)
 411 @
 412 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 413 grid <- expand.grid(x=seq(0.8,1.33,length.out=20),
 414                     y=seq(0.8,1.33,length.out=20))
 415 grid$z <- to.kcal(7^(1-1/(20*(log(grid$x)-log(grid$y))^2+1)))
 416 print(wireframe(z~x*y,grid,cuts=50,
 417           drape=TRUE,
 418           scales=list(arrows=FALSE),
 419           xlab=list("Vesicle Curvature",rot=30),
 420           ylab=list("Monomer Curvature",rot=-35),
 421           zlab=list("Curvature Backward (kcal/mol)",rot=93)))
 422 rm(grid)
 423 @
 424
 425
 426 \newpage
 427 \subsubsection{Length Backwards}
 428 \begin{equation}
 429   l_b = 3.2^{\left|l_\mathrm{ves}-l_\mathrm{monomer}\right|}
 430   \label{eq:length_backward}
 431 \end{equation}
 432
 433 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 434 grid <- expand.grid(x=seq(12,24,length.out=20),
 435                     y=seq(12,24,length.out=20))
 436 grid$z <- 3.2^(abs(grid$x-grid$y))
 437 print(wireframe(z~x*y,grid,cuts=50,
 438           drape=TRUE,
 439           scales=list(arrows=FALSE),
 440           xlab=list("Average Vesicle Length",rot=30),
 441           ylab=list("Monomer Length",rot=-35),
 442           zlab=list("Length Backward",rot=93)))
 443 rm(grid)
 444 @
 445 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 446 grid <- expand.grid(x=seq(12,24,length.out=20),
 447                     y=seq(12,24,length.out=20))
 448 grid$z <- to.kcal(3.2^(abs(grid$x-grid$y)))
 449 print(wireframe(z~x*y,grid,cuts=50,
 450           drape=TRUE,
 451           scales=list(arrows=FALSE),
 452           xlab=list("Average Vesicle Length",rot=30),
 453           ylab=list("Monomer Length",rot=-35),
 454           zlab=list("Length Backward (kcal/mol)",rot=93)))
 455 rm(grid)
 456 @
 457
 458
 459 \newpage
 460 \subsubsection{Complex Formation Backward}
 461 \begin{equation}
 462   CF1_b=1.5^{CF1_\mathrm{ves} CF1_\mathrm{monomer}-\left|CF1_\mathrm{ves} CF1_\mathrm{monomer}\right|}
 463   \label{eq:complex_formation_backward}
 464 \end{equation}
 465
 466 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 467 grid <- expand.grid(x=seq(-1,3,length.out=20),
 468                     y=seq(-1,3,length.out=20))
 469 grid$z <- 3.2^(grid$x*grid$y-abs(grid$x*grid$y))
 470 print(wireframe(z~x*y,grid,cuts=50,
 471           drape=TRUE,
 472           scales=list(arrows=FALSE),
 473           xlab=list("Vesicle Complex Formation",rot=30),
 474           ylab=list("Monomer Complex Formation",rot=-35),
 475           zlab=list("Complex Formation Backward",rot=93)))
 476 rm(grid)
 477 @
 478 <<fig=TRUE,echo=FALSE,results=hide,width=7,height=7>>=
 479 grid <- expand.grid(x=seq(-1,3,length.out=20),
 480                     y=seq(-1,3,length.out=20))
 481 grid$z <- to.kcal(3.2^(grid$x*grid$y-abs(grid$x*grid$y)))
 482 print(wireframe(z~x*y,grid,cuts=50,
 483           drape=TRUE,
 484           scales=list(arrows=FALSE),
 485           xlab=list("Vesicle Complex Formation",rot=30),
 486           ylab=list("Monomer Complex Formation",rot=-35),
 487           zlab=list("Complex Formation Backward (kcal/mol)",rot=93)))
 488 rm(grid)
 489 @
 490
 491
 492
 493
 494
 495 % \bibliographystyle{plainnat}
 496 % \bibliography{references.bib}
 497
 498
 499 \end{document}