by: David Ouyang and Dr. Jonathan Silberg

(You may also see the full spread of this article in a PDF.)

Abstract

Recent advancements in molecular biology and biochemistry allow for a new field of bioengineering known as synthetic biology. Using biological parts discovered in the last thirty years and mathematical models grounded in physical principles, synthetic biology seeks to create biological systems with user-defined behaviors. The major focus of research in this emerging field is the characterization of genetic regulation and the abstraction of biological systems to clearly defined logic circuits. With the abstraction of individual DNA sequences to known biological functions, synthetic biologists seek to create a standard list of interchangeable biological parts as the foundation of this emerging field. Through genetic manipulation, these parts are expected to be useful for programming biological machines that process information, synthesize chemicals, and fabricate complex biomaterials that improve our quality of life.

Genomic Era and Tools of the Trade

On June 26, 2000, President Bill Clinton and Prime Minister Tony Blair, along with Francis Collins, director of the Human Genome Project at the NIH, and Craig Venter, president of Celera Genomics, announced the arrival of the genomic era with the sequencing of the first draft sequence of the human genome. With this wealth of information, scientists and policy-makers alike were eager to welcome in the genomic era of genetics. Doctors dreamed of personalized medicine, where genomic information can be used to diagnose individual predispositions to cancer and disease. Politicians pondered the implications of genetic profiling, where insurance companies can potentially use genetic information to screen policyholders. The genomic era is bright with promise and unprecedented potential but also rife with social implications and practical applications.

While a sequenced genome provides a boon of new information and the scientific community is quick to emphasize the potential of this plethora of information, there are still many challenges in its interpretation and analysis. The interpretation of genomic data requires both high throughput techniques, such as microarray analysis, and heuristic algorithms in bioinformatics to analyze large amounts of data. Microarray analysis allows researchers to understand differential expression of many different proteins between different species, ages, and diseases states. With more than four billion base-pairs in the human genome and over thirty thousand open reading frames, the sheer size of the human genome requires the use of ad hoc analytical methods. The status quo approach in analyzing individual enzymes and molecules is complemented by a recent desire to understand entire systems, regulatory networks, and gene families. Exponentially increasing information on biological organisms and increasing computational power has broadened the perspective of current biological research.

Although genomic sequences provide insight into the enzymes that make up an organism, understanding of how these parts work together to produce complex phenotypes is the focus of current research. Understanding the regulation of gene expression and multicellular development will require a deeper analysis of how transcription and stability of mRNA is regulated in response to the environmental stimuli. Despite the age old debate between nature vs. nurture, it is the interplay of the environment and gene products that determine disease states and merge to create the fascinating output of life. Greater understanding of the regulation of gene products is required in determining their effects on physiology and development. Synthetic biology seeks to understand and apply understanding of biological regulation to tackle general problems.

Recombinant DNA technology laid the foundation for manipulation of biological systems on a molecular level, but recent advances in DNA sequencing and synthesis technology have greatly expanded the potential of biological engineering projects. The decreasing cost of oligonucleotide synthesis as well as improved techniques of combining oligonucleotides allows unparalleled flexibility in synthesizing long DNA sequences. From traditional methods of subcloning using restriction endonucleases and ligases to polymerase-based techniques such as gene Splicing by Overlap Extension (gene SOEing), researchers have unprecedented power in their ability to alter and characterize DNA. We can now identify new genes or regulatory sequences in diverse systems and recombine them into novel networks that attempt to recreate our understanding of existing biological systems. The rapidly expanding molecular biologist’s toolkit broadens the scope of manipulation to whole genetic systems instead of individual genes.

The current state of molecular biology has improved our understanding of the networks of biomolecular interactions that give rise to complex phenotypes and allows for unprecedented control of biological systems through clear characterization and synthetic techniques. Just as electrical engineering required increased aptness in manipulating individual circuits and transistors, biology is on the cusp of synthetic potential as new technologies overcome technical difficulties challenging previous generations of scientists.

Concept

Synthetic biology can be described as a hierarchy of fundamental biological concepts. From discrete genetic parts to whole biological circuits, each level of regulation builds upon a lower level of biological function for the ultimate goal of using biological systems to perform novel tasks or improving upon natural functions. Individual genetic parts, or particular DNA sequences with known functionality, can be integrated into genetic circuits. Genetic circuits, or new combinations of regulatory and coding sequences, can be created to produce unique behavior. Ultimately, these genetic circuits can be incorporated into biological organisms or systems.

Ongoing efforts in synthetic biology are focused on the creation of reusable, modular fragments with clearly characterized behavior and functionality in biological systems. With the discovery of the lac operon, biologists recognized the possibility for digital, discrete outputs within biological systems. Detecting the presence of lactose, the LacI repressor recognizes and binds to particular DNA sequences upstream of coding regions, regulating the transcription of the gene products in an all-or-none fashion. With clearly characterized behavior, the LacI repressor is already widely used in biotechnology applications, such as PET vectors, as integral parts of simple genetic circuits. As the biological analog of electronic circuits, researchers hope to use a growing repertoire of genetic parts to mimic logic functionalities and produce complex output.

The basic premise of synthetic biology is the ability to characterize and categorize a database of biological parts. A prominent example of this concept is the Registry of Standard Biological Parts (http://partsregistry.org/Main_Page) maintained by the BioBrick Foundation. Drawing upon the analogy of Lego bricks, synthetic biologists hope to use a standardized list of biological parts, ‘BioBricks’, to build large constructs with novel activity and unprecedented functionality. With defined activities for each component and a standardized subcloning method for combining DNA sequences, synthetic biologists hope to easily integrate ‘BioBricks’ to create novel biological circuits in a process analogous to the way computer scientists program computers. A database of DNA sequences, the Registry details the specific activity of individual sequences, the original sources and additional information necessary for synthetic biologists to integrate biological parts into particular biological systems.

Genetic Parts

There are three distinct levels where biological information can be regulated. In biological systems, information moves from DNA to RNA to protein. First proposed by Francis Crick in 1958, this “central dogma of molecular biology” addresses the detailed residue-by-residue transfer of sequential information. Synthetic biology utilizes regulatory elements at each level of this basic concept to create novel biological machines. On the DNA level, current understanding of genetic regulation reveals a complex system of promoters and terminators regulating transcription. The Registry contains a wide collection of parts for regulating transcription and translation, such as constitutive, inducible and repressible promoters, operator sequences, and ribosomal binding sites. Promoters, the 5’ upstream DNA sequencing before coding regions, determine the amount, duration, and timing of the translation. In the Registry, there is a large catalogue of terminator sequences. The 3’ region after coding regions, which form hairpin loops at the end of mRNA transcripts, cause RNA polymerase to dissociate from the template strand and end transcription. This compilation of a wide body of knowledge and literature about genetic regulation chronicles the behavior of many DNA sequences found in native systems.

Knowledge of regulation on the RNA level is applied to synthetic biology and builds upon a deep understanding of the regulation of protein production through mRNA stability and translation efficiency. Native systems display a wide range of RNA regulation that help modulate where and when particular proteins are translated. Transcribed RNA has the unique characteristic of being able to form diverse forms of secondary structure. Hairpin looping, which is intramolecular basepairing of palindromic sequences of RNA, is the basis of RNA secondary structure and can be used to create complex three dimensional structures. With the additional complexity of secondary structures, engineered RNAs function in RNA interference, as riboregulators, and catalyze key reactions. These RNA structures have been shown to mediate ligand binding and show temperature dependent activity. With temperature mediated stability, RNA sequences with designed hairpin loops can function as biological thermometers, regulating temperature sensitive expression.

The use of riboregulators is a prominent example of applying understanding of RNA behavior to regulate the expression of gene products in biological systems. Collins et al designed a system of RNA molecules that requires cooperative function of multiple RNA molecules for translation to occur. An mRNA transcript has an additional 5’ sequence complementary to the ribosomal binding site, prevent binding to the ribosome from binding to and translating the gene product. This ‘lock’ sequence can be unlocked by the regulated production of another mRNA molecule with similar homology and tighter binding affinity; allowing translation. Riboregulators with ligand mediated activity can bind to specific mRNA transcripts, silencing translation of particular genes as the result of exogenous stimuli. Both as sensors of environmental stimuli and in mediating translation, RNA has a distinct regulatory role allowing for programmable cellular behavior.

Genetic Circuits

In synthetic biology, identified regulatory components are recombined into novel networks that behave in predictable ways. An early example of a genetic circuit is the AND gate. Mimicking the functionality of digital logic of AND gates in which two unique inputs must combine to produce a positive output, Arkin and coworkers designed and modeled a genetic part to synthesize a marker protein in the presence of both salicylate and arabinose. Salicylate and arabinose are two naturally occuring, freely diffusible metabolites that bacteria normally react to; this proof of principle construct showed the ability to produce a novel reaction to simultaneous induction of both metabolites. Using two inducible promoters (NahR induced by salicylate and AraC induced by arabinose), this particular genetic part transcribed a unique T7 polymerase and the SupD amber suppressor terminator. The SupD tRNA allows translational read through at the amber stop codon, while the mutant T7 polymerase transcript includes two internal amber codons. Without the transcription of the SupD tRNA, the mutant T7 polymerase transcript would only create a nonfunctional protein product, while the SupD itself cannot induce transcription after the T7 promoter. With the combination of both gene products, a functional T7 polymerase can be expressed, which will synthesize any gene products behind the T7 promoter.

An ultimate goal of genetic manipulation is the creation of unique genetic devices or systems that can display unique characteristics or output not found in natural systems. An example of such a biological device is the repressilator, a biological device emulating the functionality of a digital oscillator which oscillates in its production of three different protein products. A system of inter-regulating gene products, the repressilator allows for sequential expression of three individual elements. Mimicking time dependent processes commonly found natural organisms, such as the KaiABC system and the circadian rhythm in photosynthetic organisms, this genetic circuit indicates the ability of simple DNA sequences to produce complex behaviors. Although this proof-of-concept constructare not as robust as natural systems, this biological device demonstrates the potential of deliberate genetic engineering to create novel output and emulate natural organisms.

A LacI repressible promoter regulates a tn10 transposon gene product which can repress another tn10 transposon promoter. This pTet promoter regulates the cI gene. A regulatory unit originally found in lambda phage, the cI protein regulates a lambda promoter that natively regulates switching between lytic and lysogenic phages in the lambda phage lifecycle. In a time dependent manner, the repressilator mimics the circadian clock found in most eukaryotic and many prokaryotic organisms.

Applications and Conclusions

In the last one hundred years, electrical systems have changed the face of the earth. Since the invention of the transistors, computers, phones, and other electronic systems have encroached upon all aspects of daily life. One can barely go through one day without use of e-mails, televisions, or cameras. Synthetic biologists dream of another world-changing revolution. Through modular parts and deliberate design, synthetic biology hopes to design biological systems to tackle challenging problems. From smart, self-regulating treatments for cancer to new solutions to the global energy crisis, the ability to engineer biological organisms has the potential to address many status quo questions. The vast natural diversity of life is a testament to the potential and opportunities available in synthetic biology.

Many different native biological organisms, such as E. coli and S. cerevesiae, are already used in many pharmaceutical and biotechnology applications. With a goal of standardization and optimization, synthetic biology allows for novel possibilities as well as improvement upon existing engineered systems. Regulating the interaction of bacteria, bacteriophage, and mammalian cells can allow for applications in medical diagnosis and treatment. The feasibility of using bacteria in biofabrication and energy generation requires designed logic functions in biological systems and biological computation. One interesting area of investigation is the removal of non-essential genes from the genome E. coli to produce an idealized minimal cell. With less chance of interfering regulatory sequences and gene products, such a minimal “cell chassis” could be the optimal shuttles for synthetic gene networks. A simplification of the cellular environment allows for greater ease in characterizing and modeling biomolecular interactions.

Utilizing the modularity of many biological systems, researchers hope to eventually produce complex behaviors through the simple combination of different biological parts. However, important considerations and research into the modularity of biological parts must still be made. In an idealized world, biological parts and coding regions could work equally well in all different cell types and organisms. Unfortunately, due to the inherent complexity of cells and intrinsically noisy nature of molecular systems, different modules might not work in different cellular environments or might not be optimized for maximum efficacy. The stochastic nature of biochemical interactions requires more work to build synthetic models and thereby understand both natural biological systems.

As genomic sequencing costs continue to decrease, the number of characterized native biological parts and unique designed parts will increase exponentially. Ultimately, synthetic biology introduces novel biological architectures not present in nature. As synthetic biology seeks to stretch the boundary of biological limits and go beyond what currently exists, questions of ethics and morality need to be addressed. What should be the limitations of investigation in this powerful field? With projects like the Venter Synthetic Genome Project, will the threshold between aggregates of molecules and life be more blurred? Should there be manipulation of the human genome, both for medicinal treatments as well as non-life threatening situations? How will intellectual property be handled, as the objects in question are inherent in natural systems? This author does not have the answers to these difficult questions, but feels that one needs to balance the potential benefits with the putative risks in this potent area of research. With great power comes great responsibility; a critical and diligent eye must be maintained in this area of active research. In addition to advancing current knowledge, it is the responsibility of the scientific community to educate the public to the potential and the risks of synthetic biology. Both to prevent Luddite reactions and to address legitimate concerns, dialogue and education are required of a field that seeks to make broad impacts on society at large.

Applying the tools and understanding of molecular biology and biochemistry, synthetic biology focuses on using current molecular tools to engineer unique biological parts and systems. Through such an engineering approach, synthetic biology also seeks to augment current approaches toward understanding regulation. Designed structures and sequences, not unique to natural systems, can be used to understand the finer details of regulation down to the very last nucleotide. As we continue to increase our knowledge of both prokaryotic and eukaryotic regulation, synthetic biologists continue to increase their repertoire of biological parts. Synthetic genome projects are currently underway and new applications such as biological computation, biological chemical fabrication, and disease treatments are being unveiled. Coupled with selection and refinement of genetic devices, deliberate genetic engineering has the potential to tackle many challenges in the near future.

David Ouyang is a sophomore Biochemistry & Cell Biology major at Baker College.

Acknowledgements

I would like to thank Dr. Jonathan Silberg for introducing me to this fascinating field of research and Dr. Daniel Wagner for his keen eye and constructive feedback. Their advice, encouragement, and support greatly aided in the writing process. I would also like to thank Erol Bakkalbasi for editing my paper.

References

1. PCR-mediated recombination and mutagenesis. SOEing together tailor-made genes. Mol Biotechnol. 1995 Apr;3(2):93-9. http://www.ncbi.nlm.nih.gov/ pubmed/7620981
2. Stephan C Schuster . Next-generation sequencing transforms today’s biology. Nature Methods – 5, 16 – 18 (2008) Published online: 19 December 2007; | doi:10.1038/ nmeth1156. http://www.nature.com/nmeth/journal/ v5/n1/full/nmeth1156.html
3. F. Narberhaus, T. Waldminghaus & S. Chowdhury. RNA thermometers. FEMS Microbiol Rev, 30(1):3-16, 2006. PMID:16438677
4. Farren J Isaacs, Daniel J Dwyer & James J Collins. RNA synthetic biology. Nature Biotechnology 24, 545 – 554 (2006) Published online: 5 May 2006; | doi:10.1038/nbt1208
5. Suess, B., Fink, B., Berens, C., Stentz, R. & Hillen, W. A theophylline responsive riboswitch based on helix slipping controls gene expression in vivo. Nucleic Acids. Res. 32, 1610–1614 (2004).
6. James J Collins et al. Engineered riboregulators enable post-transcriptional control of gene expression. Nature Biotechnology 22, 841 – 847 (2004) Published online: 20 June 2004; | doi:10.1038/nbt986
7. Poinar, H.N. et al. Science 311, 392–394 (2006).
8. Applications of next-generation sequencing technologies in functional genomics. Morozova O, Marra MA. Genomics. 2008 Nov;92(5):255-64. Epub 2008 Aug 24.
9. Molecular basis for temperature sensing by an RNA thermometer. The EMBO Journal (2006) 25, 2487–2497, doi:10.1038/sj.emboj.7601128 Published online 18 May 2006
10. A Synthetic Oscillatory Network of Transcriptional Regulators; Michael B. Elowitz and Stanislas Leibler; Nature. 2000 Jan 20;403(6767):335-8.
11. Synthetic biology: new engineering rules for an emerging discipline. Molecular Systems Biology 2 Article number: 2006.0028 doi:10.1038/msb4100073 Published online: 16 May 2006
12. Voigt C. Genetic parts to program bacteria. Current Opinion in Biotechnology Volume 17, Issue 5, October 2006, Pages 548-557 Systems biology / Tissue and cell engineering
13. “Guttmacher, Alan E., Collins, Francis S. Welcome to the Genomic Era. N Engl J Med 2003 349: 996-998
14. J Christopher Anderson, Christopher A Voigt, and Adam P Arkin. Environmental signal integration by a modular AND gate. Mol Syst Biol. 2007; 3: 133. Published online 2007 August

by: Orianna DeMasi and Kathy Li

For the full article complete with figures, please see the pdf of this article from the magazine.

Abstract

Occurrences of Lyme Disease have dramatically increased since the disease was named in the 1970’s. Much research now focuses on controlling Lyme through ticks (the vector for the disease), or White Tailed Deer (the host for ticks). Deer have recently reached surprisingly high numbers in the Northeast and it is thought that reducing deer populations will effectively decrease tick populations and thus the threat of Lyme. Consequently many towns have considered implementing ambitious yet controversial deer culling programs. As an alternative, we look at the potential for coyotes to biologically control deer populations through predation.

Coyotes, who prey on deer, have recently migrated into the Northeast from the Plains and may have been attracted to their new territory by the abundant prey supply. It is questioned whether the coyotes will act to replace the wolf as a natural control on deer population, thereby reducing Lyme Disease.

We construct a deterministic model to represent the current deer and coyote population dynamics and use this model to investigate the long term interaction of coyotes and deer. We explore the potential for coyotes to either solely act as a biological control on the deer populations or aid deer culling programs. The model is explored analytically and numerically and predicts that significant human intervention is needed to successfully control deer. Thus numerical simulations of the model and possible culling programs are provided to help highlight the system dynamics and guide culling policies.

Introduction

Recently the northeastern region of the United States has suffered an explosion of both White-Tailed Deer and cases of Lyme Disease. These two explosions are not considered to be independent and both issues greatly concern residents and policy makers of the area.

It is thought the deer explosion is a result of an increase in human population density. As humans developed the area, they created more disturbed landscapes that favor the growth of grasses, a major deer food source. Human development also conflicted with large carnivores and resulted in humans effectively killing wolves and natural ungulate predators in the region. Without predation, deer populations burgeoned in the last century [5]. Populations have gotten so large that deer have been reported in some areas at population densities as high as ten times what is thought to be healthy deer densities (20 deer per square mile vs. 200 deer per square mile) [5, 22]. Having such dense deer populations is dangerous; it leads to increased deer-car collisions [11, 24], over-grazing that can destroy natural forests, and endangerment of indigenous species [15, 4]. Possibly the greatest danger of a dense deer population, (though), is that it increases the risk of contracting Lyme disease [11, 26, 25, 22].
Lyme Disease is caused by at least three spirochetal bacteria in the genus Borrealia, but most commonly by Borrealia burgforferi. The spirochete is transmitted to individuals by ticks of the genus Ixodes. When a tick bites a host to feed on their blood, the tick transmits the spirochete to the host through its saliva and thus infects the host. Because adult Ixodes ticks frequently feed on the blood of White Tailed Deer, it is thought that more deer mean more tick and thus more Lyme. Municipalities have thus begun to explore the possibility of protecting their citizens from Lyme disease and dangers of deer overpopulation with deer control programs [5]. These programs usually propose periodically culling a portion of the deer population. Such programs have many drawbacks.

Culling is expensive and diffiult to practice as annual surveys are needed to count and monitor deer populations. Further difficulties arise with residents as the majority of land is privately owned and not all citizens morally agree with killing deer or bringing the dangers of hunting close to their homes.

We want to look for another way to potentially control deer without having to cull seasonally. We look for a biological control in the area that could act to control deer in adequate time. Specifically, we investigate if the effect of the emergent Eastern Coyote would be enough to act as this biological control.

In addition to allowing explosive deer growth, the extinction of wolves in the region has opened a niche for larger carnivores. As the emergence of a resident coyote population has occurred simultaneously, some suggest the coyote, which preys on ungulates, may be filling the wolf niche [14, 2]. Coyotes indigenous to the Great Plains and Southwest United States began migrating east and have successfully established populations in areas as far as New Brunswick and Nova Scotia [17, 2]. The eastern migration of the coyote has happened rather quickly and is thought to have resulted in a population of coyotes that is drastically different from the western coyote. The eastern coyote is physically larger [23, 13, 18, 8, 12] and has a larger home range than its western counterpart [27]. It is unknown why coyotes have moved so quickly to the Northeast, but it is believed that coyotes were attracted by abundant prey [23].

Very little is known about New England’s newest resident carnivore. It is necessary to learn about coyotes in the Northeast, as their new habitats of suburbia and New England deciduous forest differ greatly from their original home of expansive plains or the arid Southwest and puts coyotes in close proximity to humans. However, mammals with large home ranges are incredibly diffcult to track and study, so other methods of study are needed to investigate and understand coyotes on a broad scale in the Northeast [9]. These studies are especially needed if coyotes present a possible predator for deer and municipalities are looking for a method to control deer. Currently there seems to be no literature that looks at potential impact from the growth of the coyote population on Lyme. Further it appears no research has looked at exploiting coyotes as a new natural predator and control on deer populations.

We hope that employing a mathematical model might be a good way to study coyote dynamics and overcome the problems of tracking coyotes. We hope to gain insight into whether coyotes will have a significant impact on bringing deer down to acceptable densities and if the coyotes can achieve this feat within some reasonable time frame. As humans have a low tolerance of coyotes, we are also interested in if the coyotes can accomplish deer reduction with low coyote densities. Since we find mathematically that coyotes alone will not successfully control deer, we want to investigate the different types, impacts, and efficiencies of culling programs that towns might pursue. We hope to find a program that will allow us to suggest a minimally invasive cull to minimize expenses by only culling the minimum number of deer.

Modeling

Preliminary Assumptions and Notations
For simplicity, we ignore spatial variation and focus on a one square mile spatially homogeneous closed system. While this prevents modeling the varied distribution of fauna across the landscape, it will allow us to make estimates for larger scales such as multi-state domains. We also use similar assumptions to average predation and growth over a year. There is naturally some annual variability in growth and predation as birth rates and deer vulnerability vary on seasonal conditions such as snow depth [16, 17]. However this assumption is necessary to maintain an autonomous system of equations for analysis; it will also not inhibit our goal of studying the long term population dynamics.

We wish to model the density of deer D(t) and coyote C(t) populations with respect to time t in months. We begin with a classical Lotka-Volterra predator prey system,

where r is deer growth rate, a is the proportion of deer \that die in coyote-deer interactions, e is the energy that coyotes get from each killed deer, and d is the death rate for coyotes. Note this model is too simplistic, as in the absence of coyote predation, this model assumes exponential growth for deer and, in the absence of deer prey, exponential decay for coyotes. Therefore we introduce terms to more realistically portray coyote and deer interaction and reliance on the environment.

Predation
As in the classical model, we assume that at reasonable densities, deer die only as a result of coyote-deer interactions. In our sample area, this assumption is nearly reasonable as other natural predators of deer (wolves and large cats) are virtually extinct in the area. To represent this density-dependent predation, we use a Holling type II functional response term [10].

Classical mass-action predation terms such as aDC show that as prey increase, the number of prey killed by each predator increases. This is accurate if prey populations are relatively small. If prey become dense, mass action says that each coyote will kill proportionally large numbers of deer. It is more realistic that each coyote will have a maximum number of deer it can eat each month, regardless of how gigantic a population may be. Whether each coyote’s predation is maximized depends on if there are ample deer.

Holling type terms allow us to model such a cap for individual coyote predation (see Figure 1). A Holling type III term is frequently used to model mammal predation as it shows prey switching and lower predation rates if primary prey abundance is low [1]. However as studies carried out by Patterson (2000) indicate, the coyote-deer interactions better resembles a type II response

where α is the maximum deer that a single coyote will consume in one month and β controls how quickly the predation reaches α.

Note that a problem with Holling predation is that it does assume deer must reach an infinite population before individual coyotes maximize their prey consumption. Finally, we modify the Holling term to take into account the density of predators, to include C in coyote-deer interaction term.

Prey Growth
To model deer growth, we choose to use a logistic term:

where Rd is a natural growth rate for the deer population and Kd is a carrying capacity for deer. (The issue of how to chose the carrying capacity is extremely difficult, important to both coyotes and deer, and discussed later.) While the deer populations have exploded in a relatively short time, it does not seem reasonable to assume that deer continually reproduce exponentially. Furthermore, if we were to assume exponential growth and run our model with initial values as large as current population estimates, then the D(t) would quickly overrun the model and inhibit any sort of reasonable study. Recent field surveys also show that deer populations seem to be stabilizing, though only at extremely high densities [7].

Predator Growth
We want the coyote population to be correlated with deer and to grow with increased deer, which the Holling term allows. However if we use the Holling term with the extremely high populations of deer that have been observed in the Northeast, then the model will predict exploding coyote population. As there are reports of human-coyote encounters and coyotes killing pets, coyotes are often seen as a threat by humans and it is reasonable to assume that humans will control coyotes [3, 27]. If the coyote population gets too large, humans will begin killing coyotes. It is reasonable to assume that the density of coyotes that humans will tolerate is below the density of the coyotes that the current prey population could support. To impose this anthropogenic growth barrier, we multiply predation by a logistic growth term with carrying capacity Kc. The importance of Kc will become apparent later.

We also want to consider that coyotes prey on a variety of forest creatures and fruit: they are not solely dependent on the presence of deer [21, 16]. Hence we introduced a term for coyote growth independent of deer:

This term allows Rc to represent the growth from coyotes preying on species other than deer. Multiplying by a logistic term (1 − Kc ) represents the human barrier imposed on this term too. Then in the absence of deer, the coyote population will grow extremely slowly.

Carrying Capacity
The issue of carrying capacity is a delicate issue that arises quite frequently in biological problems. Mathematically, the carrying capacity represents a value which, if the density exceeds this value, growth becomes negative and the population decays back to the carrying capacity. Biologically, carrying capacity is the maximum density of a species that a given habitat can support long-term in ideal conditions.
Mathematically, carrying capacities are typically constants because they approximate long term dynamics. However, many factors affect biological carrying capacity such as food supply, climate, and over crowding. Carrying capacities can vary with respect to what is considered an external factor, such as human development, or an internal factor, such as overabundant inhabitants overgrazing and damaging the environment, thus diminishing their food source. Due to different factors, carrying capacities are an extremely delicate matter; we take Kd, Kc to be constant within our model, however, and find reasonable values from literature that predict maximum population levels.

Specialized Coyote-Deer Model
Taking into consideration all of the terms developed above, we end up with the following model:

To aid our analysis, we nondiminsionalize the system with the following substitutions:

These substitutions give us a system that behaves in the same manner and has the same equilibria as (1) but has fewer visible constants and simpler computations. The nondimensionalized system is

The reader must be cautious and recall while looking at our results that we will work with a dimensionless system. Results are dimensionless and should be analyzed as such.

Analysis

Equilibrium
We are interested in the long term interaction of coyotes and deer, so we look for points of equilibrium. More specifically, we are interested in whether the coyotes can control the deer, i.e. the existence of an equilibrium below the deer carrying capacity that has neither coyotes nor deer extinct. Equilibria occur at the intersections of the nullclines

There are six intersections and thus equilibria for our system. Given as (x, y) they are

Mathematically, coyotes controlling deer translates to an equilibrium with 0 < x< 1 and 0 <y ≤ 1. As all parameters are positive we immediately see the last equilibrium (8) with negative x value is biologically meaningless and we ignore it. It turns out that (7) does satisfy our conditions. However, to show this we must analyze the stability of the remaining relevant equilibria.

Stability
To show stability of equilibria (3) -(7), we evaluate the Jacobian at each equi-librium and examine the signs of eigenvalues. The 2 × 2 jacobian matrix is

At (3) and (5) the eigenvalues are positive and therefore equilibrium (3) and (5) (mutual extinction and coyote extinction, respectively) are unstable.
Evaluating the Jacobian at (4), we obtain

As the eigenvalues of the above matrix are negative, depending upon parameter values, we see that this equilibrium is conditionally unstable. As (4) represents coyotes at carrying capacity and deer extinct, we assume this equilibrium is unstable. Thus we force one negative and one positive eigenvalue and we obtain the parameter restriction

It can easily be shown that this restriction forces equilibrium (6) and (7) to have positive and negative x-values, respectively. We discard equilibrium (7), as biologically meaningless and concentrate on equilibrium (6), the only equilibrium with positive values for both x and y. We expect this equilibrium to be stable. Evaluating the Jacobian yields an upper triangular matrix

Therefore, J11 and J22 are the eigenvalues and must both be negative for the equilibrium to be stable.

where

Squaring the both sides of the above equation we obtain

where the rst inequality uses the restriction 9.

The above inequality is equivalent to

The three coecients of delta in the numerator of J11 can be written as

Hence,

Consider the six terms in the numerator of J11 that are independent of delta. Further, consider that due to energy transfer through tropic levels, ecosystems support more herbivores than carnivores; thus, it is reasonable that Kc < Kd.With this assumption two of the terms in the numerator are

The remaining four terms can be written as

By noting that and combining (11) and (12), we obtain

Therefore, J11 is negative and one of the eigenvalues is negative.

Now we check the sign of the other eigenvalue

Again as parameters are positive, the denominator is positive. Due to the negative sign in front, it only remains to show the sum of terms in the numerator of J22 is positive. Since there is only one negative term in the numerator, we simply show that combined with the other positive terms, the result is positive. Looking at the last three terms in the numerator and using (10) from above, we have

As the last three terms in the numerator of J22 are nonnegative, the entire numerator is positive and J22 is negative.

In conclusion, both eigenvalues J11 and J22 are negative and equilibrium (6) is stable. This means that our coyote-deer system has a stable equilibrium below the maximum number of deer and thus predicts that coyotes will have some controlling effect. Now we numerically examine the extent of this effect.

Numerical Investigation

Setting Parameters
We consider ranges of the parameters which have been gathered from the liter-ature and refine these ranges by choosing values which best model the historical growth of the two populations. The values

gave us results which best correlated with the historical dynamics described: deer population suffered until the 1950’s but exploded by the 1990’s [22] and coyote observations were sparse in the late 1950’s but more common around the 1970’s [6, 19]. This can be seen in Figure (2).

Model Results
We start the model with deer and coyotes at their respective carrying capacities, as these are the presumed current levels, and use MapleTMsoftware to numeri¬cally solve and plot our system. We consider evolution of the system over the next fifty years in Figure 3. We see that without predation, deer stay at their carrying capacity but with coyote predation, deer population drop to the level of (6), the previously found equilibrium. This shows us that coyotes do and will have an effect on deer populations.

The thin horizontal line in Figure 3 represents 20 deer per square mile which is a third of the current deer population levels and is thought to be a natural, healthy level that deer existed at in the presence of the wolf and that deer can continue to exist at without overpopulating the area [15, 4]. This is a level that is recommended by state officials for deer control programs [22, 11]. To begin to control Lyme Disease, municipalities would even like to see deer below this recommended level [5].

Deer Carrying Capactity Kd
As previously discussed, we take carrying capacity to be constant. In Figure 3 Kd is set to current level of deer populations, 60 deer per square mile. However, we can also set Kd to what the literature suggests is a healthy environmen¬tal carrying capacity, 20 deer per square mile. At this level, the logistic term causes the deer population to naturally fall to the healthy environmental carrying capacity of 20 deer per square mile within 300 months without any human intervention or predation. Even with an artificially high initial condition of 10 times the natural carrying capacity or 200 deer per square mile and no coyote predation, the deer fall to reasonable densities fairly quickly as see in Figure 4. This does not reflect reality as deer populations have remained well above the healthy carrying capacity of 20 deer per square mile and have not fallen nor are they showing signs of steep decline. Thus clearly, we cannot use the theoretical carrying capacity but must use the current population levels and let Kd=60.

Additional Coyotes
It is evident from Figure 3 that at their current densities of .2 coyotes per square mile, coyotes will not be able to control deer population to the desired levels. If we are correct in assuming that humans control coyote population size to a certain density Kc, then we can consider what would happen if humans would allow more coyotes to inhabit the area. Perhaps, if coyotes were fostered in the area to act as biological control on deer, we would see reduced numbers of deer.

We study this scenario by raising the value of Kc from .2 to .38, a doubling of the population from current levels in Massachusetts [27]. Figure 5 shows that even doubling the coyote population, while having a larger effect on the deer, will not effectively get deer to or maintain deer at desired levels.

We see that on their own, and at higher levels, coyotes will not control deer population. Even with larger populations (which are unrealistic due to the dense population in the area who are already reluctant to share the landscape with the carnivore), coyotes cannot control the deer. Hence, we now consider other efforts that would have to be taken to control deer and add the effect of culling to our model.

Culling

Analytically, we found that coyotes would indeed control the deer population to some extent. The level that coyotes could keep the deer at are below the observed densities of 60 deer per square mile as illustrated in Figure 4. However, within the range of realistic parameters, our model does not show that coyotes are capable of controlling deer to or below the 10 deer per square mile densities that are presumed optimal for controlling Lyme disease or reducing the nuisance of deer. Thus, we establish that external forces, such as active culling, will be needed to control deer and we look at some of the types of culling.

We begin our study by reconsidering deer growth and modifying our original equation to take into account a loss of deer due to an external active cull. We consider two types of culling: proportional and lump.

We consider the culling strategy that takes a proportion of the population by introducing μ, the proportion of deer killed. Frequently this sort of culling is used for modeling, as it allows for simpler analysis than the more realistic constant value or lump culling.

Continuous Proportional Culling
We consider the culling strategy that takes a proportion of the population by introducing μ, the proportion of deer killed. Frequently this sort of culling is used for modeling, as it allows for simpler analysis than the more realistic constant value or lump culling. The dimensionless culling model is

Figure 6 shows the sensitivity of the model to coyote predation. Killing only 1.2% of the deer population has a drastic effect of reducing deer. However this model also illuminates the importance of coyotes. With the increased pressure of culling, deer are more susceptible to the effect of coyote predation.

However this model seems unrealistic as culling 1.2% of deer is an extremely small portion of the deer. It is projected that much larger culls will be needed to control deer populations. We propose that the fallacy of this model is due to it assuming that the culling is continuous. This means the model is continuously taking 1.2% of the deer which is clearly unrealistic. Culling programs tend to be administered through seasonal culls or taking deer at only one period in the year and not continuously. We would like to take this seasonal culling into consideration, and thus we need a way to make the effect of culling on our model discrete.

Discrete (Seasonal) Culling
To study a realistic seasonal, proportional culling model, we utilize numerical approximating packages in MapleTM . We write a simple program which runs the system for a year, stops and subtracts a proportion μ of the deer population, and then uses that value as an initial condition for the next year. This program outputs a plot of projected yearly deer population values.
We see that with discrete culling much larger proportions of deer need to be killed in order to control populations (Figure 6). This is the more realistic result that we hoped to achieve with our MapleTMprogram. However one problem remains: time. Looking at the horizontal axis we see that culling a set proportion takes too long to get deer to or below desired levels.

Conclusion

We looked at the potential of the eastern coyote as a natural biological control on deer populations. Analytically, we found one interesting stable equilibrium in which the coyotes control the deer. Numerically analyzing this equilibrium for reasonable parameters, we find that although the coyote do have some effect on controlling deer populations, it is not enough to keep deer populations below the desired levels (20 deer per square mile). Therefore, we investigated potential anthropogenic control through culling strategies. We found that modeling continuous culling was unrealistic so we wrote a program to represent seasonal culling. This annual culling model gave far more realistic predictions yet still showed that it would take too long to control the deer if the proportion culled each year was kept constant. In conclusion, we proposed the following plan: Cull deer at high rates for a short time period, say five years, and then at more modest rates to maintain healthy populations. The predictions for such a strategy can be see in Figure 7.

Acknowledgements

The authors acknowledge with gratitude support for this work provided by the NSF REU Site Grant at Texas A&M University DMS-0552610 as well as their advisors Dr. Jay R. Walton and Dr. Yuliya Gorb.

Orianna DeMasi is a senior majoring in mathematics at McGill University. Kathy Li is a sophomore majoring in Mathematical Economic Analysis at Brown College.

References

1. Allen, L. ”An Introduction to Mathematical Biology”, 2006.
2. Ballard, W., H. Whitlaw, S. Young, R. Jenkins, and G. Forbes Predation and survival of white-tailed deer fawns in north central New Brunswick, J. of Wildlife Management, Vol. 63, pp 574-579, 1999.
3. Berchielli, L.,Impacts of Urban Coyotes on People and Pets in New York State, Proceedings of Wildlife Damage Managment Conference, 2007.
4. DeCalesta, D. Effect of white-tailed deer on songbirds within managed forests in Pennsylvania, J. Wildlife Management, Vol. 58, pp 711-718, 1994.
5. Deer Alliance. Retrieved July, 15, 2008. Web site: http://www.deeralliance.ie/
6. Fener, H., J. Ginsberg, E. Sanderson, M. Gompper Chronology of Range Expansion of the Coyote, Canis latrans, in New York, Canadian Field-Naturalist, Vol. 119, No. 1, pp 1-5, 2005.
7. Gregonis, M. 2006/2007 Aerial deer survey indicates stable population, Connecticut Wildlife, Vol. 27, No. 3 pp3, 2007.
8. Gompper, M. The ecology of northeast coyotes: current knowledge and priorities for future research, Wildlife Conservation Society Working Paper, Vol.17, pp 1-47, 2002.
9. Gompper, M., R. Kays, J. Ray, S. Lapoint, and D. Bogan, A Comparison of Noninvasive Techniques to Survey Carnivore COmmunities in Northeastern North America,Wildlife Society Bulletin, 2006.
10. Holling, C.S. Components of predation as revealed by a study of small-mammal predation of the European Pine Sawfly, Canadian Entomologist, Vol. 91, pp 293-320, 1959.
11. Kilpatrick, H. J., and LaBonte, A. M. ”Managing Urban Deer in Connecticut: a guide for residents and communities”, Hartford: Bureau of Natural Resources, 2007.
12. Kyle, C. J., Johnson, A.R., Patterson, B.R., Wilson, P.J., Shami, K., Grewal, S.K., and White, B.N. Genetic nature of eastern wolves: past, present, and future, Conservation Genetics, Vol. 7, pp 273-287, 2006.
13. Larviere, S. and Crete, M. The size of eastern coyotes (Canis latrans): a comment, J. of Mammalogy, Vol. 74, pp 1072-1074, 1993.
14. Mathews, N. E., and Porter, W.F. Maternal defense behavior in white-tailed deer, pp 141-160 in A. H. Boer, editor. ”Ecology and manage¬ment of the Eastern Coyote. University of New Brunswick”, Fredericton, New Brunswick, Canada, 1992.
15. Mc Shea, W.J., and Rappole, J.H., Managing the abundance and diversity of breeding bird populations through manipulation of deer populations, Conservation Biology, Vol. 14, No.4, pp 1161-1170, 2000.
16. Patterson, B.R., Benjamin, L.K., and Messier, F. Prey switching and feeding habits of eastern coyotes in relation to snowshoe hare and white-tailed, Canadian Journal of Zoology, Vol.76, pp 1885-1897, 1998.
17. Patterson, B.R. and Messier, F. Factors influencing killing rates of white-tailed deer by coyotes in eastern Canada, J. Wildlife Management, Vol. 64, No. 3, pp 721-732, 2000.
18. Peterson, R.O., and Thurber, J.M. The size of eastern coyotes ( Canis latrans): a rebuttal, J. of Mammalogy, Vol. 74, pp 1072, 1993.
19. Pringle, L.P. Notes on coyotes in Southern New England, J. Mammal, Vol.41, pp 278, 1960.
20. Rand, P.W., Lubelczyk, C., Holman, M.S., Lacombe, E.H., and Smith Jr., R.P. Abundance of Ixodes scapularis (acari: Ixodidae) after complete removal of deer from an isolated offshore island endemic for Lyme disease, J. Medical Entomology, Vol. 41, pp 779-784, 2004.
21. Samson, C., M. Crete,Summer food habits and population density of coyotes, Canis latrans in boreal forests of southeastern Quebec, Canadian Field-Naturalist, Vol. 111, No. 2, pp 227-233, 1997.
22. Stafford, K.C. 2004. ”Tick management handbook: an integrated guide for homeowners, pest control operators, and public health officials for the prevention of tick-associated disease”, The Connecticut Agricultural Ex-periment Station, New Haven, Connecticut, 2004.
23. Thurber, J. M. and Peterson,R. O. Changes in body size associated with range expansion in the coyote ( Canis latrans), J. of Mammalogy, Vol. 72, pp 750-755, 1991.
24. Walter, W.D., Perkins, P.J., Rutberg, A.T., and Kilpatrick, H.J. Evaluation of immunocontraception in a free-ranging suburban white-tailed deer herd, Wildlife Society Bulletin, Vol. 30, pp 186-192, 2002.
25. Wilson, M.L., Ducey, A.M., Litwin, T.S., Gavin, T.A., and Spiel-man, A. Microgeographic distribution of immature Ixodes dammini ticks correlated with deer, Medical and Veterinary Entomology, Vol. 4, pp 151¬159, 1990.
26. Wilson, M.L., Adler, G.H., and Spielman, A. Correlation between abundance of deer and that of the deer tick, Ixodes dammini (Acari: Ixodi¬dae), Annuals of the Entomological Society of America, Vol. 7, pp 172-176, 1985.
27. Way, J.W., Ortega, I.M., and Auger, P.J. Eastern coyote home range, territoriality, and sociality on urbanized Cape Code, Northeast Wildlife, Vol. 57, 2002.

by: Arindam Sarkar and Dr. John S. Olson

For the full article complete with figures, please see the pdf of this article from the magazine.

Abstract

Last summer, Dr. Mitchell Weiss and his colleagues at Children’s Hospital of Philadelphia discovered a new hemoglobinopathy in a baby from Toms River, New Jersey, who was born cyanotic and with enlarged spleen and liver tissues. Sequencing of the baby’s hemoglobin alleles revealed a missense mutation in a segment of DNA that codes for the gamma chains of fetal hemoglobin (HbF), the oxygen-carrying protein in red blood cells of human fetuses. The objective of our work is to use recombinant DNA technology to construct the Hb Toms River mutation, γ Valine 67 (E11) Methionine, in plasmid DNA which can then be used to express and purify mutant protein using E. coli. We plan to characterize the mutant HbF in order to understand its clinical manifestations and, perhaps, to develop treatments options. This paper provides an overview of HbF developmental biology, our initial hypothesis of how the Hb Toms River mutation might lead to cyanosis, and our strategy for expressing and characterizing the γ Val67 to Met mutation in recombinant HbF.

Introduction

During a consultation with pediatricians at Children’s Hospital of Philadelphia, Dr. Mitchell Weiss discovered a new blood disorder in a child who was born cyanotic and with an enlarged spleen and liver. These symptoms resolved roughly two months after her birth, and she was normal and healthy by six months. Based on their initial clinical observations, Dr. Weiss and the treating physicians suspected that a mutant fetal hemoglobin might be the cause of the baby’s symptoms, so they drew small amounts blood samples for analysis from the baby several days after birth. They discovered that her condition appears to fall into a class of hematological disorders known as hemoglobinopathies, which are genetic defects in the DNA sequences that produce hemoglobin. Hemoglobin is the primary oxygen transport protein in humans, and when the baby’s DNA was analyzed, it was discovered that a Val67 to Met (V67M) mutation was present in one of the child’s γ chain alleles. This mutation occurs in a region of DNA that gives rise to the eleventh amino acid along the E helix of the γ globin chain, which is called Val (E11) for its spatial location in the three dimensional structure of hemoglobin subunits (Figure 1).

The original mutant fetal protein could not be studied directly because of physical and clinical limitations which prevent the withdrawal of significant amounts blood from an anemic infant. Another problem was that fetal hemoglobin production switches to adult hemoglobin production shortly after birth as part of normal developmental processes. Thus, resolution of the cyanotic condition occurred when the γ gene (characteristic of HbF) was silenced, and only normal adult hemoglobin was present in the baby’s red blood cells, which occurred 6 to 8 weeks after birth. To obtain enough starting materials for study, we chose to produce mutant HbF in our laboratory using recombinant technology. The objective of our work is to use structural biology to characterize the γ V67M mutation in HbF, examine the role of the E11 position in O2 binding in γ chains, and then understand why the mutation caused cyanosis and spleen enlargement.

Hemoglobin Development

Hemoglobin is a complex iron-containing protein in the blood that picks up oxygen from the lungs and carries it to respiring cells; at the same time, it assists in transporting carbon dioxide away from the peripheral tissues. Mammalian red cell hemoglobins are tetramers consisting of four polypeptide chains and four planar prosthetic groups known as heme molecules [2, 3, 9]. Each red blood cell contains about 280 million hemoglobin molecules [7].

Different kinds of hemoglobins are commonly identified by the specific combination of polypeptide chains or subunits within each tetramer. During development and birth, three main types of hemoglobin are expressed (Figure 2). The first type is known as embryonic hemoglobin, which consists of two α and two γ globin chains. The low oxygen conditions of the uterine wall demand a higher oxygen affinity than either HbF or adult hemoglobin confer, but embryonic hemoglobin functions well in this environment. After 10 to 12 weeks of development, the primary form of hemoglobin switches from embryonic hemoglobin to fetal hemoglobin HbF (α2γ2). At this point, the fetus’s red blood cells have access to the oxygen passing though the placenta and umbilical cord. Like embryonic hemoglobin, HbF has a higher oxygen affinity than adult hemoglobin, thus allowing the fetus to extract oxygen from the mother’s blood in the placenta.

When a baby is born and begins to breathe air, γ chain production ceases and β chains are produced which result in adult hemoglobin HbA (α2β2) production. At birth, HbF comprises 50-95% of the child’s hemoglobin. These levels decline to almost zero after six months as adult hemoglobin synthesis is completely activated. Hemoglobin genes exist on Chromosomes 11 and 16 (Figure 3).

Genetic abnormalities can suppress the switch to adult hemoglobin synthesis, resulting in a condition known as hereditary persistence of fetal hemoglobin [6]. In adults, HbF production can be rekindled pharmacologically, which is one of the main treatment options for sickle-cell disease [5]. The mechanisms by which erythroid cells switched from the synthesis of HbF to that of HbA during the neonatal period appeared normal for the patient with HbF Toms River. As a result, this genetic disorder did not persist as a threat to the child.

Analysis of the V67M Mutation

A necessary stage for understanding the Hb Toms River disorder will be acquiring large amounts of the cyanotic child’s mutated fetal hemoglobin. In this case, our only choice was to construct the mutation in vitro with recombinant DNA techniques and then express the mutant γ chain with wild-type α chains in bacteria. The intestinal bacterium, Escherichia coli, is an excellent choice for expressing recombinant proteins because of its high tolerance for synthesizing large amounts of heterologous proteins and the ease of performing site-directed mutagenesis on plasmids that can be taken up by this bacterium. The plasmid system pHE2 was originally developed by Chien Ho’s group at Carnegie Mellon University (Shen et al., 1993) to produce adult hemoglobin, and we obtained the pHE2 expression system for HbF from Professor Kazuhiko Adachi at Children’s Hospital of Pennsylvania (Adachi, 2002). This vector contains one wild-type α gene, along with one wild-type γG gene from human HbF. We created the single-site V67M mutation using the Stratagene QuikChange Site-Directed Mutagenesis Kit. The plasmids for expression of hemoglobins were transformed into E. coli JM109 cells. E. coli cells were grown in 2x YT medium. Expression was induced by adding isopropyl-β-thiogalactopyranoside (IPTG) at 0.1 mM at 37oC and then supplemented with hemin (30 μg/ml). The harvested cell lysate was passed through a Zn2+ binding column, Fast Flow Q-Sepharose column, and finally a Fast Flow S-Sepharose column using an FPLC.

We are currently evaluating the purity and authenticity of our HbF mutant by performing gel electrophoresis and protein sequencing reactions from aliquots of the purified protein. Our long-term goal is to characterize HbF Toms River in terms of its relative stability and O2 affinity with the hope that recombinant technology can help us understand the clinical symptoms of the hemoglobinopathy and perhaps suggest a treatment.

For the past twenty years Dr. John Olson’s laboratory in the Department of Biochemistry & Cell biology at Rice has been examining O2 binding to mutants of mammalian myoglobin and the α and β subunits of HbA. Much of this work has focused on amino acid substitutions within the oxygen binding pocket, including at the valine E11 position. In 1995, an undergraduate honors research student in Dr. Olson’s laboratory, Joshua Warren (Rice BA 1996; Yale PhD 2002), used sperm whale myoglobin (Mb) as a model system to examine the effects of valine E11 to methionine, phenylalanine, tyrosine, and tryptophan mutations on oxygen binding. All four of these amino acids have much larger side chains which fill up the interior portion of the pocket which captures diatomic gases, including O2, CO, and NO. Warren observed dramatic decreases in the rates of O2 uptake and release due to the valine E11 to methionine replacement. Similar marked decreases compared to wild-type Mb were observed for the Phe, Tyr, and Trp mutations, which also decrease the size of the binding pocket [4].

The mechanism of O2 binding to either Mb or a Hb subunit is analogous to catching a baseball in a fielder’s glove. As the thumb opens, by upward and outward movement of the histidine E7 side chain (Figure 4), incoming oxygen can be “caught” in the pocket of the glove. If the available space of the glove is made too small by limiting it with a large amino acid like methionine, the ball or O2 will “bounce” back out of the globin requiring that multiple tries be made until it is finally captured and bound to the iron atom. Thus, we expect the valine E11 to methionine mutation to appreciably slow O2 binding.

At the moment, a structure of the γ valine E11 to methionine mutant has only been simulated (Figure 4) and recombinant HbF Toms River not been characterized. However, based on Warren’s results with Mb, we predict that O2 binding may be slowed so much that red blood cells containing the HbF mutant cannot uptake oxygen quickly enough during passage through the placenta or new born lungs. Consequently, the blood will only be partially saturated with O2 and appear a purplish or cyan color associated with cyanosis. We are currently setting up the HbF expression system and, as a control, show that wild-type γ subunits have kinetic and stability parameters very similar to those of HbA β subunits.

With luck this system will allow complete characterization of the HbF Toms River mutation, and recombinant DNA technology will allow us to study the clinical disorder associated with this hemoglobinopathy without having to obtain the infant’s blood. This in vitro approach represents an important advance in characterizing genetic defects in hemoglobins and will provide a general approach for determining the underlying mechanisms for the phenotype associated with the hemoglobin mutations and possible treatments to restore normal physiological function.

Arindam Sarkar is a sophomore double-majoring in Biochemistry & Cell Biology and Policy Studies at Lovett College.

Acknowledgements

I thank Dr. John S. Olson and Todd Mollan for their useful comments, Ivan Birukou for a helpful discussion on the purification of recombinant hemoglobin, Dr. Jayashree Soman for her expedient provision of images, and once again Todd Mollan for his tireless support and tutelage.

References

1. Adachi, K., Zhao, Y., and Surrey, J. (2002) Assembly of Human Hemoglobin (Hb) beta – and gamma -Globin Chains Expressed in a Cell-free System with alpha -Globin Chains to Form Hb A and Hb F. J. Biol. Chem. 277, 13415-20.
2. Bunn, H.F. and Forget,B.G. (1986) Hemoglobin: Molecular, Genetic and Clinical Aspects. Saunders, Philadelphia, PA.
3. Dickerson, R. E., and Geis, I. (1983) Hemoglobin: Structure, Function, Evolution, and Pathology, pp. 21-26, Benjamin -Cummings Publishing Co., Menlo Park, CA
4. Nienhaus, K., Deng, P., Olson, J.S., Warren, J.J., and Nienhaus, G.U. (2003). Structural Dynamics of Myoglobin: Ligand migration and binding in valine 68 mutants. J. Biol. Chem. 278, 42532-42544.
5. OS Platt. (2008) Hydroxyurea for the treatment of sickle cell anemia. N. Engl. J. Med. 358, 1362–9.
6. S Friedman, E. Schwartz. (1976) Hereditary persistence of foetal haemoglobin with beta-chain synthesis in cis position. Nature. 259, 5539.
7. Sears, Duane W. 1999. Overview of Hemoglobin’s Structure/Function Relationships. February 2005.
8. Shen,T.J., Ho, N.T., Simploaceanu, V., Zhou, M., and Ho, C. (1993) Production of unmodified human adult hemoglobin in Escherichia coli. Proceedings of the National Academy of Science 90, 8108-8112.
9. Steinberg, M.H., Forget, B.G., Higgs, D.R., and Nagel, R.L. Disorders of Hemoglobin . Cambridge University Press, Cambridge, 2001.

by: Sarah Wu and Professor Yousif Shamoo

Abstract

Only recently has the U.S. public become aware of the dangers of antibiotic resistance. Overuse of antibiotics is chiefly responsible for the proliferation of super-resistant strains. An intimate understanding of resistance mechanisms will prove crucial for continued success against pathogenic bacteria. Traditionally, antibiotics have targeted processes such as cell wall synthesis and DNA replication. However, the decreasing effectiveness of optimizing existing antibiotics has prompted scientists to explore new strategies, such as using combinations of antibiotics, examining evolutionary pathways leading to resistance, and targeting metabolic processes.

Introduction

The problem of antibiotic resistance catapulted to the attention of the U.S. public when seventeen-year-old Ashton Bonds died in October 2007 from MRSA (methicillin-resistant Staphylococcus aureus), which he contracted from his high school locker room.¹ Studies conducted by the Centers for Disease Control and Prevention estimate deaths in the U.S. due to MRSA in 2005 exceeded those of HIV-AIDS, Parkinson’s disease, emphysema and homicide.² While traditionally associated with hospital ICUs and nursing homes, deadly bacterial infections emerged in schools, gyms, and day cares, alarming the nation.

Historically, antibiotics have demonstrated a remarkable effectiveness against bacteria. Many consider antibiotics as one of the most important discoveries in modern science. The discovery of penicillin by Alexander Fleming in 1928 prevented many deaths and amputations on the battlefields during WWII.³ Sulfa drugs not only accelerated research in the pharmaceutical industry but also were successful starting points for the creation of new drugs to treat a variety of diseases.⁴ Despite their amazing effectiveness, antibiotics are slowly losing their edge against pathogens, as seen by the recent rise of strains such as MRSA, which doubled from 127,000 infections in 1999 to 278,000 in 2005 and increased from 11,000 to more than 17,000 deaths.⁵

As medicine integrated antibiotics into common practice as a way to treat a variety of infectious diseases such as bronchitis, pneumonia, and syphilis, antibiotics have become the expected and preferred cure in the vocabulary of the general public for unpleasant symptoms. These drugs have come to be seen as the magic cure-all that doctors can easily dispense to their patients, who are increasingly expectant of instant treatments. According to rationalmedicine.org, physicians may over-prescribe antibiotics due to fear of lawsuits of negligence if drugs are not prescribed. In addition to facing pressure from drug companies to prescribe their medications, physicians risk losing their patients to other physicians who are less hesitant to prescribe antibiotics.⁶ Sore throat is usually caused by a viral infection (and is therefore unresponsive to antibiotics) with only 5-17% of cases caused by bacteria. However, a 2001 study by Linder and Stafford revealed that in the period from 1989 to 1999, there was a decrease in the use of recommended antibiotics like penicillin for treating sore throat and an increase in the use of nonrecommended broad spectrum antibiotics, such as extended-spectrum macrolides and fluroroquinolones.⁷ Another study showed that half of the antibiotic prescriptions written by emergency medicine physicians were for viral infections.⁸ Because of a variety of factors leading to faulty administration, antibiotics are not being used effectively.

Antibiotics are also widely used in agriculture. Frequently animals that are raised for consumption in the U.S. come from factory farms that maximize yield by placing as many animals as possible into a limited space. Close quarters such as these facilitate the spread of infectious disease, making it necessary to treat the animals regularly with antibiotics as a preventative measure. Antibiotics are also administered to make up for poor conditions like inadequate nutrition and imperfect caretaking.⁹ This is concerning as there have been reports of animal to human transmission of disease. A 2006 study found a case where a strain of MRSA from a family’s pig farm had infected a baby girl, demonstrating the ease with which resistant strains can arise in farm animal populations and spread to humans.¹⁰

Resistance against antibiotics is not a novel phenomenon in biology and has in fact occurred for millions of years. Streptomyces, a genus of bacteria that lives in soil and water, excretes a variety of anti-microbial substances into the nearby soil to prevent the growth of competing microorganisms.¹¹ Penicillin mold has similar defenses. It is not surprising that we derive many of our antibiotics from nature.

If bacteria are quickly becoming resistant to antibiotics, why do pharmaceutical companies hesitate to develop replacement drugs? In short, drug development is very difficult, expensive, and unprofitable. It is estimated that the cost of developing one successful drug is around $802 million. Drugs take about 12 years to get approved by the FDA and only about one out of 5,000 substances are approved, a 0.02% success rate. The pharmaceutical company then has the rights to the drug patent for 20 years before the drug becomes generic, allowing anyone to produce it.¹² In addition, the scope of bacteria against which most antibiotics are effective quickly narrows with increased usage. The drug company is under pressure to start selling as much of the drug as possible to recoup the cost of development, which means aggressive marketing to physicians. Physicians are subsequently faced with a dilemma, as they want to limit the use of antibiotics to maintain their effectiveness yet feel pressured into prescribing them. Currently, a large portion of research is dedicated to the mechanisms and selection of antibiotic resistance, which scientists attribute to evolutionary processes.

Evolution’s role in antibiotic resistance

Evolution is the continual adaptation of a population to its environment. Mutations at the nucleotide level influence the structure of the proteins to produce a phenotype that could confer increased resistance to antibiotics. When a naive population of bacteria is exposed to an antibiotic, the vast majority is eliminated. However, the surviving bacteria are resistant and will give rise to antibiotic-resistant progeny. Thus, repeated administration of the same antibiotic results in the proliferation of resistant phenotypes and decrease in antibiotic efficacy. Because bacteria are unicellular organisms that divide roughly every thirty minutes, they quickly adapt to environmental challenges such as antibiotics. The combined effect of many divisions and mutations gives rise to numerous resistance mechanisms.

Antibiotics: Mechanism and Response

Bacterial cell walls serve the critical function of maintaining osmotic stability, making their synthesis a favorite target of bactericidal antibiotics. Bacteria are classified into two categories based on the type of cell walls that they have: gram-positive or gram-negative (Fig. 1). Gram refers to the staining protocol developed by Hans Christian Gram in 1884.¹³ The crystal violet stain used in this procedure will color the outer peptidoglycan wall of gram-positive bacteria purple but leaves the gram-negative bacteria pink, as the peptidoglycan net is contained inside its dual membranes. This peptidoglycan net is made up of alternating units of N-acetylmuramic acid (NAM) and N-acetylglucosamine (NAG) which are cross-linked by the enzyme transpeptidase, also known as penicillin-binding proteins (PBP). The peptidoglycan net is a critical component to the bacteria as it maintains the high osmotic pressure inside the bacteria.

β-lactam drugs (Fig. 2), such as penicillin and ampicillin, are the most widely used antibiotics. β-lactams are a class of drugs that inhibit bacterial cell wall synthesis. Because their structure mimics the penultimate D-Ala-D-Ala pentapeptide chain of NAM, they function by being mistakenly used by the transpeptidases as a substrate, resulting in the acylation of the enzyme. The acylated transpeptidase is unable to hydrolyze the β-lactam and is thus not functional, which hinders cell wall synthesis. The cell wall becomes susceptible to the autolytic enzymes that degrade the cell wall and eventually becomes permeable to the environment, leading to bacterial lysis.¹⁴

Bacteria have responded in three main ways to β-lactam drugs. The active site of the transpeptidase can be changed to have a lower affinity for β-lactams. In gram-negative bacteria, the cell can also restrict the influx of drugs into the periplasm by altering the expression of outer membrane proteins to restrict drug entry or by the expression of antibiotic efflux pumps such as MexA.¹⁵ The most common mechanism of β-lactam resistance uses an enzyme to inactivate the drug before it can bind transpeptidase. This enzyme, known as β-lactamase (Fig. 3), is very similar in structure to the PBP and uses a strategically placed water molecule to hydrolyze the β-lactam ring of the drug, rendering it inactive. However, as new forms of β-lactam drugs have been introduced, bacteria have evolved more efficient forms of this enzyme.

Cell wall synthesis, however, is not the only common target of antibiotics. Fluoroquinolones (Fig. 2) inhibit DNA replication by stabilizing the enzyme-DNA complex created by the enzymes topoisomerase II (DNA gyrase) or topoisomerase IV, blocking DNA synthesis. Bacterial mechanisms for resistance to fluoroquinolones are very similar to those for β-lactam resistance. Mutations to the target sites of DNA gyrase and topoisomerase IV decrease drug binding, increasing expression of efflux pumps to remove the drug. Alternatively, losing porins in the membrane restrict entry of the drug into the cytoplasm.¹⁶

Current Research in the Field

For many years, the main response to developed resistance was to introduce a new form of the antibiotic. Unfortunately, little progress has been made recently in developing new effective compounds, eliminating the possibility of a quick score in terms of discovering a bacterial growth inhibitor.¹⁷ This necessitates the study of novel approaches to eradicating pathogenic bacteria.

Drug cocktails

Some researchers have pursued the idea of combining antibiotics in hope of a synergistic effect. The combination of two drugs is expected to lead to a more effective therapy, as is the case with the cocktail of drugs given to HIV patients. However, predicting the efficacy of drug combinations has proven to be difficult, making direct experimentation necessary for verifying efficacy. Cefotaxime and minocycline were found to be an effective combination against Vibrio vulnificus¹⁸, although the treatment of sepsis with the addition of aminoglycosides to β-lactams was found to have no advantage and actually increased risk of nephrotoxicity.¹⁹ One study characterized combinations of drugs as synergistic, additive, or antagonistic (depending on if the cumulative effect of the drugs is greater than, equal to, or less than the effect of their individual activity) and demonstrated how certain combinations of drugs could select against drug resistance. By conducting a competition assay between doxycycline-resistant and wild type E. coli at certain concentrations of doxycycline and ciprofloxacin, the wild type strain was found to have a growth advantage over the doxycycline-resistant strain.²⁰ While this study was limited to sublethal concentrations of antibiotics, it suggests the potential to forestall drug resistance through combinations that select against the development of resistance.

Predicting the evolutionary trajectories of antibiotic resistance

Other research concentrates on the evolutionary processes of resistance. One study identified five key β-lactamase mutations that led to a 100,000-fold increase in resistance against cefotaxime. While in principle there are 120 (5!) possible pathways (also known as trajectories) to reach the final state of accumulating all five mutations, the study demonstrated that a large number of these pathways were already inaccessible as a beneficial mutation has a much higher probability of fixation in the population as opposed to a deleterious or neutral one. Other combinations of those 5 mutations did not increase resistance unless certain mutations preceded others sequentially. In the end, they found only 18 probable trajectories leading towards the super effective enzyme with five mutations.²¹ This study shows the surprisingly limited nature of evolutionary pathways and the possibly predictable nature of evolution. While this study was strictly theoretical in its execution, similar investigations into the evolutionary trajectories of antibiotic resistance at the molecular level are being conducted using in vivo models with the hope of understanding the mutational landscape of the development of drug resistance.

Targeting Bacterial Metabolism

Other studies have looked beyond conventional drug targets and in bacterial metabolism, suggesting that all antibiotics kill bacteria in the same way by causing them to produce hydroxyl radicals, which damage proteins, membrane lipids, and DNA. A methodical sequence of experiments probes this idea, starting with experiments that show the increased production of hydroxyl radicals in bacteria when treated with bactericidal antibiotics. The Fenton reaction, which produces hydroxyl radicals by the reduction of hydrogen peroxide by ferrous iron, is considered the most significant contributor of hydroxyl radicals. Bacteria were found to live longer when exposed to thiourea, a hydroxyl radical scavenger, as well as 2,2’-dipyridyl, an iron chelator that inhibits the Fenton reaction. To determine the source of the ferrous iron as either extracellular or intracellular, a knockout strain (∆tonB) was created with disabled iron import as well as a knockout (∆iscS) with impaired iron-sulfur cluster synthesis abilities. While ∆tonB exhibited no advantage against antibiotics, ∆iscS showed reduced hydroxyl radical formation and cell death, pointing to intracellular ferrous iron as the source of hydroxyl radicals. It is established that ferrous iron is released when superoxide damages the iron-sulfur clusters and the most superoxide formation comes from electron transport chain oxidation and conversion of NADH to NAD+. Gene expression microarrays revealed that upon exposure to antibiotics, NADH dehydrogenase I was a key upregulated pathway. Bacteria respond to hydroxyl radicals by activating RecA, which stimulates SOS response genes to initiate DNA repair mechanisms. RecA knockout strains had a significant increase in cell death compared to wild type. The study proposes the following mechanism: antibiotics stimulate oxidation of NADH, which induces hyperactivation of the electron transport chain, thus stimulating superoxide formation that damages iron-sulfur clusters in the cell. These clusters release ferrous iron, which is oxidized by the Fenton reaction, producing the hydroxyl radicals. Bacteria have SOS response mechanisms that are activated by the RecA gene, leading to the stimulation of SOS response genes. By knocking out this gene, they found the bacteria had a significant increase in sensitivity to antibiotics. The study demonstrated the possibilities of targeting the TCA cycle or respiratory chain when developing new antibiotics.²²

Complementary strategies

While a large burden rests in the hands of the scientists and drug companies, there are steps the public can take to decrease the rate of proliferation of resistant strains. The federal Center for Disease Control and Prevention recommends not pressuring one’s health care provider to prescribe antibiotics and to follow directions when taking medication by only taking antibiotics for bacterial infections as well as completing the prescribed dose of antibiotics completely to avoid sparing bacteria that could potentially become drug resistant.²³ Improving infection control would prevent the spread of resistant strains as rapid methods of identifying pathogens would lead to faster isolation of colonized patients, making it harder for disease to spread.²⁴ The simple practice of washing hands can limit the spread of disease. Citizens can also advocate legislation that limits the use of antibiotics in the agricultural industry.
As observed in the laboratory and daily life, carelessness with the application of antibiotics results in the disappearance of our painstakingly acquired advantage against bacteria. But even with cautious usage, antibiotics will still select for resistant strains. Continuous ongoing research is critical to discovering successful novel treatment strategies.

References

1. Urbina, Ian. “Schools in Several States Report Staph Infections, and Deaths Raise the Alarm.” The New York Times. 19 Oct 2007.
2. Sack, Kevin. “Deadly Bacteria Found to Be More Common.” The New York Times. 17 Oct 2007.
3. “penicillin.” Encyclopædia Britannica. 2008. Encyclopædia Britannica Online. 23 Mar. 2008 .
4. Lesch, John E. The First Miracle Drugs: How the Sulfa Drugs Transformed Medicine. New York: Oxford University Press, 2007.
5. Klein E, Smith DL, Laxminarayan R (2007). “Hospitalizations and Deaths Caused by Methicillin-Resistant Staphylococcus aureus, United States, 1999–2005”. Emerg Infect Dis 13 (12): 1840–6.
6. Kakkilaya, B.S.. “Antimicrobials: The Rise and Fall”. rational medicine.org. 03/24/08 .
7. Jeffrey A. Linder and Randall S. Stafford
Antibiotic Treatment of Adults With Sore Throat by Community Primary Care Physicians: A National Survey, 1989-1999
JAMA, Sep 2001; 286: 1181 – 1186.
8. Diane M. Birnbaumer. Emergency physicians overprescribe antibiotics for viral respiratory infections. Journal Watch Emergency Medicine, April 2006.
9. “feed.” Encyclopædia Britannica. 2008. Encyclopædia Britannica Online. 23 Mar. 2008 .
10. Huijsdens et al. Community-acquired MRSA and pig-farming. Ann Clin Microbiol Antimicrob (2006) vol. 5 pp. 26
11. “Streptomyces.” Encyclopædia Britannica. 2008. Encyclopædia Britannica Online. 24 Mar. 2008 .
12. Crowner, Robert. “Drug Pricing vs. Drug Development.” Acton Commentary. 15 May 2002.
13. “Gram stain.” Encyclopædia Britannica. 2008. Encyclopædia Britannica Online. 23 Mar. 2008 .
14. Babic et al. What’s new in antibiotic resistance? Focus on beta-lactamases. Drug Resist Updat (2006) vol. 9 (3) pp. 142-56
15. Wilke et al. Beta-lactam antibiotic resistance: a current structural perspective. Curr Opin Microbiol (2005) vol. 8 (5) pp. 525-33
16. Chen et al. Molecular mechanisms of fluoroquinolone resistance. J Microbiol Immunol Infect (2003) vol. 36 (1) pp. 1-9
17. Projan. New (and not so new) antibacterial targets – from where and when will the novel drugs come?. Current opinion in pharmacology (2002) vol. 2 (5) pp. 513-22
18. Chiang et al. Synergistic antimicrobial effect of cefotaxime and minocycline on proinflammatory cytokine levels in a murine model of Vibrio vulnificus infection. J Microbiol Immunol Infect (2007) vol. 40 (2) pp. 123-33
19. Paul et al. Beta lactam monotherapy versus beta lactam-aminoglycoside combination therapy for sepsis in immunocompetent patients: systematic review and meta-analysis of randomised trials. BMJ (2004) vol. 328 (7441) pp. 668
20. Chait et al. Antibiotic interactions that select against resistance. Nature (2007) vol. 446 (7136) pp. 668-71
21. Weinreich et al. Darwinian evolution can follow only very few mutational paths to fitter proteins. Science (2006) vol. 312 (5770) pp. 111-4
22. Kohanski et al. A Common Mechanism of Cellular Death Induced by Bactericidal Antibiotics. Cell (2007) vol. 130 (5) pp. 797-810
23. “Antibiotic/Antimicrobial Resistance Prevention Tips”. Centers for Disease Control and Prevention. 03/24/08 .
24. Lowy. Antimicrobial resistance: the example of Staphylococcus aureus. J Clin Invest (2003) vol. 111 (9) pp. 1265-73

by: Yuekai Sun and Professor Jun Lou

Abstract

Molecular dynamics (MD) is a computer simulation technique that studies particle models by stimulating the time evolution of a system of interacting particles. MD is commonly employed in materials science, nanotechnology and biochemistry to study various processes at the atomic scale. The techniques of MD have also been adapted for use in fluid dynamics and astrophysics. In this article, we introduce the basic science and computational methods behind MD and explore some of its current applications.

Given a system consisting of n particles subjected to interactions described by a given force field, can the trajectory of every particle be predicted? 18^th century French astronomer and mathematician Pierre-Simon Laplace said of solving the above problem, called the n-body problem, “Given for one instant an intelligence which could comprehend all the force by which nature is animated and the respective situation of the beings who compose it—an intelligence sufficiently vast to submit these data to analysis—it would embrace in the same formula the movements of the greatest bodies of the universe and those of the lightest atom; for it, nothing would be uncertain and the future, as the past, would be present to its eyes.”

Molecular dynamics (MD) solves the n-body problem by numerically integrating the equations of motion for every particle. There are two common forms of MD, ab-initio MD and classical MD. Ab-initio MD solves the many-body problem for a system governed by the Schrödinger equation with the Born-Oppenheimer approximation (without this approximation, the Schrödinger equation cannot be solved even numerically for many systems). At every timestep, the Schrödinger equation for electrons is solved to compute the potential. The position and velocity of every nucleus are then propagated by numerically integrating the classical equations of motion. To solve the Schrödinger equation for electrons, various approximations derived with methods such as the tight-binding model or density functional theory are used. Ab-initio MD can simulate quantum mechanical phenomena, such as bonding between atoms and heat conduction in metals, but the computational cost limits the size of systems that can be studied with this model to a few thousand atoms.

Classical MD simplifies the ab-initio approach by replacing the quantum mechanical potential computations with a parameterized empirical potential function that depends only on the position of nuclei. This potential function is either fitted to experimental data or derived from the underlying quantum mechanical principles. Classical MD can be used to study systems with up to millions of atoms but cannot simulate phenomena that depend on electron behavior. We will restrict our discussion to classical MD.

The name molecular dynamics is misleading because given the almost universal nature of the many-body problem, MD can model any system that involves interaction between particles. MD is employed in materials science to study a variety of phenomena, including the behavior of materials under stress, crack propagation, and how various defects affect the strength of materials. MD is also commonly used in biochemistry to study the behavior of macromolecules. The techniques of MD have also been generalized to thermal science and astrophysics, where particle models are used to study various hydrodynamic instabilities and phase transitions in thermal science and the structure of the universe in astrophysics.²

Interaction computations

The computation of interactions between atoms is typically the most computationally expensive segment of a MD simulation. For a pair-wise potential (a potential described by a potential function that does not account for the environment of atoms), the interaction computations scale with order O(N²). Various methods have been developed that reduce the order to O(N) for short-range potentials (potentials that decay faster with respect to distance than the dimension of the system are considered short-range) and O(NlogN) for long-range potentials. The most common methods are the linked-cell and neighbor list methods for short range potentials and Ewald summation method for long-range potentials.

Both the linked-cell and neighbor-list method assume that the potential and its spatial derivative at the location of an atom can be approximated by the superposition of the potential contributions from atoms within a certain interaction cutoff distance. Mathematically, the potential function U(r_ij) can be approximated by the truncated potential function Û(r_ij).

But to compute the distance between every atom pair costs O(N²) calculations. The linked-list method decomposes the system spatially into cells of side length greater than or equal to the cutoff distance. Every atom only interacts with atoms in the same cell or an adjacent cell so only distances between atom pairs in the same and adjacent cells must be computed. Atoms are resorted into cells every 10 to 20 timesteps. This reduces the complexity of the computation to O(N) .²

The linked-list method is commonly used in parallel molecular dynamics codes because it simplifies communications between threads. Every thread is assigned a cell and is passed a list of the atoms in its cell and every adjacent cell every time atoms are resorted into cells.

Figure 1 Linked-cell method: To compute the interactions of the filled in atom, only atoms in the same (dark-shaded) and adjacent (light-shaded) cells are considered. Only atoms within the shaded disk interact with the darkened atom.

The neighbor list method creates lists of all atoms within a certain distance of every atom so every atom only interacts with atoms in its neighbor list. The neighbor list cutoff distance must be greater than the interaction cutoff distance and the difference determines the frequency at which the neighbor lists must be recreated. The neighbor list method also reduces the complexity of the computation to O(N) (2).

Figure 2 Neighbor list method: The neighbor list for the filled-in atom contains all the atoms in the light-shaded disk. Only atoms within the dark-shaded disk interact with the filled in atom.

For long-range potentials, the interaction cutoff distance must be very large for the potential at the location of an atom to be approximated by the superposition of the potential contributions from atoms within the cutoff distance. The Ewald summation method is the most common method and reduces the computational complexity to O(Nlog(N)) ³ by decomposing the system spatially into cells and assuming the interactions of an atom with every other atom in the system can be approximated by its interactions with every other atom in the same cell and their images in periodic cells tessellated along all dimensions.

Mathematically, this can be expressed as:

where N is the number of atoms in the same cell, n_x, n_y, n_z are half the number of periodic cells in the x, y, z directions respectively, Ø(r_ij) is the part of the function that does not decay with respect to distance, L is the side length of the cell and d is the dimension of the system. The Ewald summation method then decomposes the above sum into a short-range term that can be summed using techniques for short-range potentials and a long-range term that can be summed with Fourier methods. A detailed mathematical introduction to the Ewald summation method is beyond the scope of this article and we refer the reader to ³ for a detailed introduction to the Ewald summation method.

Integration methods

The solution to an n-body problem can be chaotic and it is fruitless to seek an accurate trajectory beyond several “collision steps”. The main criteria for choosing an integration method for the equations of motion in MD are therefore energy conservation, the ability to accurately reproduce thermodynamic ensemble properties and computational requirement. Interaction computations are the most computationally expensive segment of an MD simulation and an integration method requiring more than 1 interaction computation per timestep is inefficient unless it can maintain the same accuracy with a proportionally larger timestep. This criterion disqualifies the common Runge-Kutta methods.

Two commonly used families of integration methods are the Verlet algorithm and its variations and predictor-corrector methods. The Verlet algorithm can be derived from the Taylor expansion of the position variable (r):

Adding the expressions for r(t+h) and r(t-h) and solving for r(t+h):

One drawback of the Verlet algorithm is the lack of a velocity term in the expressions, making it difficult to compute the atomic velocities. Common variations of the Verlet algorithm designed to correct this drawback are the leap-frog algorithm and velocity-Verlet algorithm.

The predictor-corrector family of methods (PC methods) all consists of 3 steps. First, atomic positions and velocities are propagated according to the Taylor expansions of the position and velocity variables. Then, the interactions are computed based on the new positions and the accelerations are compared to the accelerations predicted with the Taylor expansion of the acceleration variable. Finally, the positions and velocities are “corrected” according to the difference between the computed and predicted accelerations. Although PC methods require 2 interaction computations per timestep, they can maintain the same accuracy with a timestep that is more than twice as long.³

PC methods have mostly been discarded in favor of the simpler Verlet methods because Verlet methods give better energy conservation and are easier to implement. We refer the reader interested in a derivation of the basic PC algorithm to. ¹

Constraint dynamics

Geometric constraints can be implemented in the context of the above integration methods to simulate rigid molecules. The rigid molecule approach is not appropriate for large molecules. Methods to simulate flexible molecules have been developed although we will not introduce them in this article. We refer the reader to¹ for a detailed introduction to the simulation of flexible molecules. The most common method is the SHAKE algorithm and its variations. The SHAKE algorithm first advances the atomic positions and velocities ignoring the constraints. The new positions and velocities are then corrected by introducing a constraint force, the magnitude of which is determined by the displacement of the new positions from their constrained positions.

The k^th constraint constraining the distance between two atoms can be expressed mathematically as:

where r_i and r_j are the positions of atoms i and j and d_k is the bond length between them. These constraint equations represent forces that are described by their action on the dynamics of the system and are not derived from the potential. These forces can be expressed within the equations of motion as constraint potentials:

where n is the number of constraints and λ is a constant that must be computed in order for the constraint to be satisfied. We numerically integrate the equations of motion with the Verlet algorithm.

where is the unconstrained position of atom i. r_i and r_j must obey the constraint equations σ_k:

where r_k is the constrained bond length. The value of λ_k can be solved for with standard root-finding techniques. In the SHAKE algorithm, the system of k equations is solved with iterative methods [1].

Thermostats and barostats

MD with the classical equations of motion simulates a system at constant volume and energy. Every timestep, a new state of the system, called a microstate, is generated. Because every microstate generated has the same volume and energy, MD with the classical equations of motion can be thought of as sampling from a set of microstates with a given number of atoms, volume and energy, called a microcanonical/NVE ensemble (so called because the Number of particles, Volume and Energy are constant).

It is often desirable to sample from a canonical/NVT (Number of particles, Volume and Temperature are constant) or isothermal-isobaric/NPT (Number of particles, Pressure and Temperature are constant) ensemble to maintain compatibility with experimental studies. Various methods have been developed to constrain temperature and pressure and are called thermostats and barostats respectively. Most modify the system dynamics by modifying the equations of motion. Common thermostat algorithms include the Berendsen thermostat and Nosé-Hoover thermostat and common barostats include the Parrinello-Rahman method.

The equations of motion for most thermostat algorithms take the general form:

where γ can be thought of as a coefficient that adjusts the system kinetic energy to match the desired temperature. The effect of a thermostat algorithm on system dynamics can also be described as a scaling of internal velocities (velocities relative to center of mass) at every timestep by a factor λ. Mathematically, this can be expressed as:

where is the unscaled internal velocity and h is the timestep.

Because the temperature cannot change instantaneously, λ(t,0) = 1.

If we let in (*), then the above expression agrees with (*).

The Berendsen thermostat is based on Newton’s law of cooling:

where T₀ is the heat bath temperature and τ is an empirical coupling parameter that quantifies the strength of the coupling between the system and heat bath. The finite difference quotient is:

The velocity scaling factor at every timestep is determined by the ratio between the system temperature and heat bath temperature and the coupling parameter.

Note that the coupling parameter affects the system dynamics like a damping parameter. The Berendsen thermostat, therefore, does not generate a true canonical velocity distribution because the coupling parameter damps out system temperature fluctuations.

Both the Nosé-Hoover themostat and the Parrinello-Rahman method are based on the extended system approach. In this approach, additional virtual degrees of freedom are introduced into the system that can be used to regulate system properties.

A detailed introduction to the Nosé-Hoover themostat and the Parrinello-Rahman method is beyond the scope of this article. We refer the reader ⁴ for a detailed introduction to the Nosé-Hoover thermostat and²
for a detailed introduction to the Parrinello-Rahman method.

Conclusion

According to American physicist Richard Feynman, “all things are made of atoms and everything that living things do can be understood in terms of the jiggling and wiggling of atoms”. MD is a computer simulation technique that simulates the “jiggling and wiggling of atoms” and has been hailed as a way of predicting the future by animating nature’s forces. In this article, we introduced the basic method behind MD to provide the reader with the necessary background knowledge to understand an MD simulation.

References

[1]. Rapaport, D. C.The Art of Molecular Dynamics Simulations. Cambridge : Cambridge University Press, 2004.
[2]. Griebel, Michael, Knapek, Stephan and Zumbusch, Gerhard. Numerical Simulation in Molecular Dynamics: Numerics, Algorithms, Parallelization, Applications. Berlin : Springer, 2007.
[3]. Leach, A. R.Molecular Modeling: Principles and Applications. Upper Saddle River : Prentice Hall, 2001.
[4]. Frenkel, Daan and Smit, Berend.Understanding Molecular Simulation: From Algorithms to Applications. San Diego : Academic Press, 2002.
[5]. Hünenberger, Philippe. Thermostat Algorithms for Molecular Dynamics Simulations. Advanced Computer Simulation. Berlin : Springer, 2005, pp. 105-149.

by: Timothy Kinney and Professor Bruce Johnson

Abstract

One of the most significant developments in the field of Raman spectroscopy was the discovery of Surface-Enhanced Raman scattering (SERS). Raman signals from the early days of SERS were known to be enhanced by a factor up to a million-fold, and more recently by extraordinary factors of up to 10^12-14. To understand what this means we should understand how Raman scattering works, and then apply this understanding to the Surface-Enhanced case. SERS can be thought of as at least two kinds of enhancement: an electromagnetic enhancement which is well understood and a chemical enhancement which is not completely understood. Charge-Transfer has been shown to contribute some chemical enhancement, but more mysteries remain. We are working to characterize the wavelength dependence of the Anti-Stokes to Stokes ratio, thus deriving what contributions are made by metal-adsorbate-complex transitions approximately resonant with both the laser and surface plasmon frequencies in the absence of Charge-Transfer effects. Work by this group and by others indicates that anti-Stokes to Stokes ratios can be sensitive indicators of this chemical enhancement mechanism in SERS when it occurs.

Background: Raman Scattering

A new type of light-matter interaction was experimentally verified by C.V. Raman and K.S. Krishnan in 1928, now called Raman scattering.¹ In contrast with Rayleigh scattering, it is the scattering of photons with a change in frequency. Though the effect holds for all wavelengths of light, it is easier to consider a monochromatic case. During a Raman event, an incident photon of a particular wavelength and energy interacts with the vibronic modes of a matter system. The term vibronic indicates a combination of electronic states and vibrational states. The energy of the system increases by an amount allowed by vibronic transitions and a photon having lost this exact energy is re-emitted. The molecule thus moves from a lower energy state to a higher energy state, while photons of a higher energy are absorbed and photons of a lower energy are emitted. This is called a Stokes shift, but each vibronic transition also has an Anti-Stokes shift. (See Figure 1.) In the Anti-Stokes case, the matter is already in an excited vibrational state. An incident photon interacts with the excited system, resulting in the emission of a higher energy photon and a loss of energy in the system as it returns to a lower energetic state. This happens less often and is dependent on how many molecules in the system are in an excited state, generally due to thermal interactions. The probability of any Raman event occurring is on the order of one for every million incident photons, but it is highly dependent on the polarizability of the matter, and the frequency of the incident radiation. The ratio of Anti-Stokes events to Stokes events is a useful quantity for observing resonance between vibronic transitions and incident radiation.

Since the advent of laser technology in the mid-20th-century, scientists have had access to collimated beams of nearly monochromatic radiation. From this, Raman spectroscopy was developed to probe the vibronic transitions of matter systems by examining the emitted Raman light with a spectrometer. If the energy of the incident radiation is known, the energy of the emitted radiation can be measured, and the energy gained by the system can be calculated exactly. This energy then corresponds to vibronic modes in the system, offering experimental verification of theoretical calculations. Like Infrared or Fluorescent Spectroscopy, Raman Spectroscopy identifies groups of molecules by their chemical bonds, but also provides additional information. The vibrational modes of matter are easily characterized and studied, and a large body of work dedicated to this field exists. In Raman spectroscopy, vibrational modes are called Raman modes. To display Raman spectra, the wavenumber of the incident radiation (usually a laser) is marked as zero. In what is known as the wavelength shift, the wavenumber of each Raman mode is shifted from zero, corresponding to the energy of that mode. For example, a carbon-carbon stretch occurs at approximately 1600 wavenumber shift. It will always occur at 1600 wavenumber shift, regardless of the frequency of the incident radiation. Changing the incident radiation changes the absolute wavenumber position, but not the relative shift. Note that varying the surrounding environment or the physical state of the matter may shift the Raman modes slightly. Figure 4 shows two Surface-Enhanced Raman spectra of 1-dodecanethiol. Several carbon-carbon modes are visible between 1050 and 1200 wavenumber shift.

Background: Surface-Enhanced Raman Scattering

Due to the rarity of Raman events, the Raman signal is often very weak relative to the intensity of incident radiation. The experimentalist desires very intense incident radiation to be ensured of a strong signal. However, over-heating the sample can lead to unwanted side effects, so a balance must be sought between strength and authenticity of the signal. One method to mitigate this is to place the matter in close proximity with a noble metal, such as gold, silver, or copper. Because the intensity of Raman events is dependent on the polarizability of the matter, the strength of the signal is closely dependent on the number of free electrons in the system. Noble metals have an abundance of free electrons that become polarized in resonance with the incident radiation and the molecules being studied, resulting in an increased intensity of Raman photons. This is called Surface-Enhanced Raman Scattering, or SERS. When light interacts with a noble metal, the free electrons on the surface tend to polarize in resonance with the light, and this collective motion is called a surface plasmon. The peak resonance of the surface plasmon depends on the type of noble metal and its geometry. Understanding surface plasmons and their wavelength dependence is critical to understand surface enhancement.

In the late seventies, electrochemists took spectra of pyridine adsorbed on silver, copper or gold electrodes electrolyte solutions. ^3-5 Scanning the potential of the electrode (by comparison with Saturated Calomel Electrode), they demonstrated variations in the Raman signal. It was discovered that oxidation/reduction cycling would dissolve and redeposit monolayers of metal, augmenting the SERS enhancement by two orders of magnitude. This additional enhancement was attributed to molecular roughening of the surface of the electrodes.⁶ The SERS intensity was found to increase with the number of deposited monolayers for small numbers of layers.⁷ Allen et al. were able to show that multiple mechanisms contribute specified magnitudes to overall enhancement. They suggested that approximately 3.5 x 10¹ enhancement could be contributed from a surface-roughness mechanism, and 2 x 10³ enhancement could come from roughness-independent mechanisms for pyridine adsorbed on copper electrodes.⁸

Wavelength Dependence of SERS: Multiple Mechanisms

Pockrand et al. determined that the multiple mechanisms were not limited to substrate activity. They collected Raman spectra of pyridine on evaporated silver films and on silver electrodes, and compared SERS-active films with SERS-inactive films investigating wavelength dependence of Raman intensities. ⁹ They explicitly noted that excitation profiles could be plotted showing SERS intensities dependent on the laser excitation wavelength, with a single maximum in the visible region for each vibrational mode. (See Figure 2.) Their paper demonstrated that all the excitation profiles they took have similar resonance behavior, but that the curves shift to greater energy of the incident light with the increasing energy of the vibration. The surface plasmon enhancement of silver was qualitatively the same with evaporated silver films and with silver electrodes, ruling out mechanisms solely dependent on the substrate. Regarding the adsorbed molecule, they found that doublet signals were obtained from thick layers of pyridine on SERS-active films; one of the pair correlated to the bulk pyridine and the other to a red-shifted version of the same vibrational mode in the surface pyridine. This corroborates the conclusion that the SERS response for pyridine is dependent on both surface plasmon activity attributed to the silver substrate and excitation energy of the laser. Furthermore, the large enhancement factors from SERS-active films could not be described by local field effects alone, and Pockrand concluded that a chemical contribution must play a role.

The multiple mechanism picture has prevailed over the years, leading to the proposal of various models, more or less dependent on specific data obtained from SERS experiments. In general, it is now well-accepted that a local-field electromagnetic enhancement occurs, which resonates with vibronic modes in the molecule and thus substantially increases the intensity of Raman scattered radiation. This mechanism requires that the molecule’s Raman-shifted vibronic frequencies fall within the spectral width of the substrate surface plasmons, although direct proximity between metal and molecule are not required. ^{10, 11} It is also well-accepted that this mechanism does not completely describe SERS. For physisorbs/chemisorbs, molecule-specific chemical effects are found in enhancement patterns and intensities. ¹⁰ A Charge-Transfer mechanism, requiring direct proximity between metal and molecule, has been suggested to explain these.^11-13 Charge-Transfer involves the overlapping of molecular orbitals to provide pathways for the migration of electrons from the metal to the adsorbate, or vice versa. This is thought to provide enhancement by modifying the Raman polarizability tensor, causing increased intensity in Raman scattered radiation, though some scientists argue that the Charge-Transfer mechanism should be regarded as an altered adsorbate complex rather than an enhancement mechanism.¹¹

Early work with electrodes, colloids and island films was insufficient to fully understand the complexity of SERS substrate plasmon resonance activity. Researchers called for experiments with increased control of surface roughness, but technology was lacking. ⁸ A solution to this problem has been developed at Rice University and involves the fabrication of nanoscale particles, like nanoshells. The study of nanoparticles has led to major advancements in the understanding of SERS, such as Jackson et al. recently investigating the relationship between nanoshell geometry and SERS enhancement. ¹⁴ Using nanoshells, a new monodispersed substrate, Jackson et al. showed that, in particular for those modes which the laser is resonant with, the Raman Effect is enhanced proportionately with the density of the nanoshells. ¹⁴ In other work by the same authors, it was discovered that SERS intensity was linearly dependent on the density of the nanoshell substrate, proving that the SERS response in their experiments was due to single nanoshell resonance response, not dimer resonances. ¹⁵ It is important to note that tuning the laser wavelength to excite single nanoshell resonances or dimer resonances results in data corroborating one theory or the other. This reinforces the notion that wavelength dependence is critical to understanding SERS mechanisms.

Maher et al. specifically investigated wavelength dependence of SERS on multiple substrates by plotting the anti-Stokes to Stokes ratio for p-mercaptoaniline (pMA) at different excitation energies and discovered that asymmetries arise. ¹⁶ At a fixed temperature a calculable amount of anti-Stokes emission is usually expected based on the thermal activity of the sample. They used the ratio of anti-Stokes to Stokes intensities as a sensitive indicator of resonant conditions involving both the metal and the adsorbate molecule. They argued that this asymmetry is systematically structured according to underlying resonances that need more than a single wavelength to be seen. One series of data they took for laser line 633 nm produced a higher anti-Stokes to Stokes ratio than expected, and another series using a 514 nm excitation produced a lower ratio than expected. (See Figure 3.) This may correspond to an anti-Stokes emission intensity peak around the 585 nm to 610 nm region. Their data lacked the excitation energy density necessary to determine this with certainty, but they did show that the anti-Stokes to Stokes ratio is a potentially valuable means of mapping the resonant interaction between the probes and the metal.

Advancing Our Understanding

We anticipate that in some cases unexpected resonances (called “hidden” by Maher et al. ¹⁶) for the metal-molecule-complex will be observed, even in the absence of Charge-Transfer processes. We are specifically interested in 1-dodecanethiol, though we also intend to do experiments with p-mercaptoaniline, due to the other work that has been conducted at Rice University. ^10,16 Initial observations for the 1-dodecanethiol and gold colloids system using a Renishaw Invia Raman Microscope have been made. (See Figure 4.) We are currently building a Tunable Raman Darkfield Microscope system to restrict the region of measurements to those where the gold colloids have deposited in a smooth and dense manner. Using a dye laser, we can tune continuously from 558 nm to 606 nm, providing good coverage of the visible spectrum in the region where gold colloids enhance strongly. Compiling excitation profiles as a function of excitation energy, we will calculate the theoretical thermally-allowed anti-Stokes to Stokes ratio for the sample and then examine deviations from this result. Building on the work of Gibson et al. (see subsequent section of this paper), ¹⁰ we may prove or disprove the theory of hidden resonances for the metal-molecule-complex in the absence of Charge-Transfer for pMA attached to gold nanospheres and nanoshells. This work can later be applied to other substrates and molecules to provide data to support the general case.

Theoretical Support

As theoretical support for our work, we turn to Gibson et al. who independently ran theoretical density-matrix calculations for pMA and obtained results corroborating the notion that unexpected resonant behaviors between molecule-metal-complex may play a significant role in Surface Enhancement. ₁₀ Their model was able to obtain anti-Stokes resonant excitation behavior in the absence of Charge-Transfer, supporting our hypothesis that these unexpected resonances occur when the anti-Stokes emission is in resonance with vibronic transitions of the molecule-metal-complex. Gibson et al. suggest that anti-Stokes emissions which are resonant with a metal-adsorbate molecule may provide anti-Stokes intensity peaks which then fall away in intensity on either side of the peak. Their calculations agree with the work of Maher ¹⁶ and may explain why his observed anti-Stokes to Stokes ratio is asymmetrical to thermal expectations; lower or higher, depending on laser excitation energy.

Acknowledgements

The author would like to thank Rice University and the Rice Quantum Institute. He also acknowledges funding provided by the National Science Foundation and hopes that Research Experience for Undergraduates programs will continue to be funded for the long-term development of research sciences.

References

1. Long, D.A. “Raman Spectroscopy”, (McGraw-Hill, Great Britain, 1977).
2. Image from Vrije University Brussels, 2007, Dec 18, . Accessed 2008 Mar10. Axis labels added.
3. M. Fleischmann, P.J. Hendra and A.J. McQuillan, Chem. Phys. Letters 26, (1974) 123.
4. U. Wenning, B. Pettinger, H. Wetzel, Chem. Phys. Lett. 70, (1980) 49.
5. J. Billmann, G. Kovacs and A. Otto, Surface Sci. 92, (1980) 153.
6. D.L. Jeanmarie and R.P. Van Duyne, J. Electroanal. Chem. 84, (1977) 1.
7. B. Pettinger and U. Wenning, Chem. Phys. Letters 56, (1978) 253.
8. Craig S. Allen, George C. Schatz, Richard P. Van Duyne, Chem. Phys. Letters 75, (1980) 201.
9. I. Pockrand, J. Billman, A. Otto, J. Chem. Phys. 78 (11): 6384-6390 1983.
10. J.W. Gibson, B.R. Johnson, J. Chem. Phys. 124, 064701 (2006).
11. M. Moskovits, Rev. Mod. Phys. 57 (3): 783 1985.
12. M. Osawa, N. Matsuda, K. Yoshii, I. Uchida, J. Phys. Chem. 1994, 98 12702-12707.
13. J.R. Lombardi, R.L. Birke, T. Lu, J. Xu, J. Chem. Phys. 84 (8) 4174, 1986
14. J.B. Jackson, S.L. Westcott, L.R. Hirsch, J.L. West, N.J. Halas, Appl. Phys. Lett. 82, 257 (2002).
15. J.B. Jackson, N.J. Halas, Proc. Nat. Acad. Sci. USA 101, 17 930 (2004).
16. R.C. Maher, J. Hou, L.F. Cohen, E.C. Le Ru, J.M. Hadfield, J.E. Harvey, P.G. Etchegoin, F.M. Liu, M. Green, R.J.C. Brown, M.J.T. Milton, J. Chem. Phys. 123, 084702 (2005).