1 COMPUTATIONAL INVESTIGATIONS OF ORGANIC MATERIALS FOR HYBRID NANODEVICE AND OPTOELECTRONIC APPLICATIONS By JASMINE DAVENPORT CRENSHAW A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2011
2 2011 Jasmine Davenport Crenshaw
3 To my Lord and Savior, Jesus Christ, all that I am is because of you. To my loving husband, Anthony Crenshaw; my supportive parents, Clarence & Segrid Davenport; my loving sisters, Jamila Davenport and Torrie Davenport Walker; my extended family and friends, I could not have completed this journey without you. You have graciously been my support system. I also dedicate this work to my grandparents, the late Thermon Hudson, Geneva Hudson, Clarence Davenport, and Ethel Davenport. Even though none of you graduated from high school you instilled in me the importance of education. I stand on t he shoulders of giants!
4 ACKNOWLEDG E MENTS To God be the glory for the great things He has done! Luke 12:48 states, "To whom much is given, much shall be required". I am so honored for all that He has birthed out of me through this experience. I know th at I am nothing without you. To my darling husband, Anthony Marcus Crenshaw, I want to say, "I love you"! You stuck by me through thick and thin, and provided wisdom and encouragement during times that I felt I could not go on any more. You sacrificed so m uch for me during this process, and I am forever grateful. To my parents, Clarence and Segrid Davenport, you are my solid rock and always have my back! To my sisters, Jamila Davenport and Torrie Davenport Walker, you kept me humbled and showed me the impor tance of balance with family and work. I love you both! To my mother n law, Maria Crenshaw, Uncle and Aunts, John Covington, Sandra Covington, and Marisa Casado, thank you for your contin ued love, support, and advice. I appreciate you always having my best interest at heart! To Drs. Simon Phillpot, Henry Hess, and Susan Sinnott, thank you for your inspiration, guidance, and mentorship. There were times when I wanted to throw in the towel but you were patient and believed in me. I am forever thankful. To Drs ., Laurie Gower, Richard Dickinson, and Josephine Allen, I am grateful to have this experience with you. You all have been extremely helpful and challenged me in numerous areas t o utilize my maximum potential. To the faculty and staff of North Carolina Agr icultural and Technical State University, I am appreciative of my foundation and training. To Dr. Cesar Jackson, thank you for being my role model. You encouraged me to pursue a doctoral degree. To Drs. Solomon Bililign, Claude Lamb, Rita Lamb, Floyd James Samuel Danagoulian, and James Williams you invested countless hours of teaching and mentorship. You
5 equipped me with the tools needed for research in graduate school. To the Ronald E. McNair Postbaccalaureate Achievement Program at NC A&T State Universit y, under the leadership of Dr. Jane Brown and Mrs. Geraldine Burnette, the National Institute of Health Minority Access to Research Careers Undergraduate Student Training in Academic Research (MARC U STAR), under Drs. James Williams and Claude Lamb, and th e National Science Foundation Talent 21 program, under the direction of Dr. Cesar Jackson, Mrs. Sunnie Howard sharpened my skills and introduced me to research. To Dr. Henry Frierson, Dean of the Graduate Scho ol at the University of Florida, thank you for your dedication and support to students of color. You kept me focused and personally held me accountable to obtain my doctoral degree. I appreciate your leadership and being a part of my academic "posse". To t he University of Florida Office of Graduate Minority Programs and Ronald E. McNair Program, your offices have provided encouragement and support. I personally would like to thank Dr. Laurence Alexander, Mr. Earl Wade, Mrs. Sarah Perry, Mrs. Janet Broiles, and Mrs. Edna Daniels. I truly consider your offices my adopted family. To Drs. Nedialka Iordanova and Tzvetelin Iordanov, Michelle Morton, and W. Joseph Barron from Georgia Southwestern State University, thank you all for the computational methods trainin g sessions and all of the helpful discussions. I must recognize all of the past and current research group members of Drs. Simon Phillpot and Susan Sinnott for their valuable group interactions, friendships and support. My sincere gratitude is extended to Drs. Tao Liang, Rakesh Behera, and Priyank Shukla for their invol vement in my research projects.
6 I must recognize the love and support of my family and friends: the Crenshaws; Davenports; my spiritual parents, Pas tor George and Lady Michele Dix; Covingtons ; Casados; Harrisons; Parkers; Staples; Ruffins; Thorpes ; my church, PASSAGE Family Church (Gainesville, FL) ; Mt. Gilead Baptist Church (Durham, NC) ; Dr. Charlee and Andrew Bennett; Dr. Danyell Willson; Dr. Samesha Barnes; Dr. Tara Washington; Dr. Tameka P hillips; Dr. Beverly Brooks Hinojosa; Anderson Prewitt; Wanda Eugene; Sydni Credle; Cephas Small; Tamekia Jones; Darina Palacio; William Hardy; Fatmata Barrie; Krystal Lee; Crystal Moore; Angela Henderson; Michael and Krystal Anthony ; Clyde Hicks; Ruth Law rence; the late Helen Daniel; the late Senator Jeanne Hopkins Lucas Last but not least, I thank the National Science Foundation DMR 0645023, National Science Foundation Alliance for Graduate Education in the Professoriate (SEAGEP), National Consortium for Graduate Degrees for Minorities in Engineering and Science, Inc. (GEM), the Eastman Kodak Company, and the University of Florida Ronald E. McNair Graduate Assistantship Program for providing financial assistance to fund my graduate education.
7 TABLE OF C ONTENTS page ACKNOWLEDGEMENTS ................................ ................................ ............................... 4 LIST OF TABLES ................................ ................................ ................................ .......... 10 LIST OF FIGURES ................................ ................................ ................................ ........ 11 ABSTRACT ................................ ................................ ................................ ................... 13 CHAPTER 1 GENERAL INTRODUCTION ................................ ................................ .................. 15 2 KINESIN POWERED MOLECULAR SHUTTLE ................................ ..................... 18 2.1 Introduction ................................ ................................ ................................ ....... 18 2.2 The Micr otubule Protein Filament ................................ ................................ ..... 19 2.3 Fluorescent Microscopy ................................ ................................ .................... 20 2.4 Hybrid Design Model ................................ ................................ ......................... 20 3 CE LLULAR AUTOMATON SIMULATION METHOD ................................ .............. 28 3.1 Biomolecular Shuttle Transport ................................ ................................ ......... 28 3.2 Previous Computational Models Examining Microtubule Dynamics ................. 29 3.2.1 Fluorescent Microscopy ................................ ................................ ........... 29 3.2.2 Off lattice Monte Carlo Simulation ................................ ........................... 29 3.3 Cellular Automaton ................................ ................................ ........................... 31 3.3.1 Background on Cellular Automaton ................................ ......................... 31 3.3.2 Types of CA methods ................................ ................................ .............. 32 3. 4 Fortran 90 ................................ ................................ ................................ ......... 33 3.5 Development of the Simulation Tool ................................ ................................ 34 3.5.1 Overview ................................ ................................ ................................ 34 184.108.40.206 Definition of lattice ................................ ................................ .......... 34 220.127.116.11 Simulation parameters ................................ ................................ ... 35 18.104.22.168 Description of simulated microtubules ................................ ............ 37 22.214.171.124 Motion of microtubules ................................ ................................ ... 37 126.96.36.199 Turning of microtubules ................................ ................................ .. 38 3.5.2 Binding Interactions ................................ ................................ ................. 39 3. 5.3 Treatment of Nanospool Formations ................................ ....................... 40 3.6 Simulation Memory Consumption ................................ ................................ ..... 41 3.7 CA Algorithm Verification ................................ ................................ .................. 41
8 4 MT SIMULATION RESULTS AND DISSCUSSION ................................ ................ 52 4.1 Scenario 1: Validation of the Simulation Tool ................................ ................... 53 4.1.1 Persistence Length/Turn Probability ................................ ........................ 54 4.1.2 Validation of Nanospool Formation ................................ ......................... 55 4.1.3 Number of Spools Generated and Nanospool Circumferences ............... 55 4.1.4 Nanospool Circumference Radial Distribution ................................ ...... 57 4.2 Scenario 2: Effects of Simulation Parameters on Nanospool Formations ......... 58 4.2.1 Simulation Parameters ................................ ................................ ............ 59 4.2.2 Nanospool Formation Versus Time ................................ ......................... 60 4.3 Scenario 3: Application of the Simulation Tool in the Dependence on Trajectory Persistence Length ................................ ................................ ............. 61 4.3.1 Introduc tion ................................ ................................ .............................. 61 4.3.2 Experimental Results ................................ ................................ ............... 63 4.3.3 Discussion ................................ ................................ ............................... 66 4.4 Summary ................................ ................................ ................................ .......... 70 5 POLYTHIOPHENE THIN FILM GROWTH VIA SURFACE MODIFICATIONS ....... 84 5.1 Thin Films ................................ ................................ ................................ ......... 84 5.2 Growth of Organic Films ................................ ................................ ................... 84 5.3 Organic Material Surface Modification ................................ .............................. 85 5.4 Surface Polymerization by Ion Assisted Deposition ................................ .......... 87 6 HYBRID FUNCTIONAL METHODS ................................ ................................ ........ 90 6.1 Basic Quantum Mechanics to the Schrdinger Equation ................................ .. 90 6.1.1 Time Independent Schrdinger Equation ................................ ................ 90 6.1.2 Born Oppenheimer Approximation ................................ .......................... 90 6.2 Approaches to Approximate Solutions in the Schrdinger Equation ................. 91 6.2.1 Hartree Product ................................ ................................ ....................... 91 6.2.2 Hartree Fock (HF) Approximation ................................ ............................ 92 6.2.3 Density Functional Theorem (DFT) Methods: Hohenberg Kohn and Kohn Sham ................................ ................................ ................................ ... 94 188.8.131.52 Hohenberg Kohn theorems ................................ ............................ 94 184.108.40.206 Kohn Sham equation ................................ ................................ ..... 95 6.2.4 Relative Advantages and Disadvantages of HF vs DFT methods ........... 96 6.3 Hierarchy of Approximations ................................ ................................ ............. 97 6.3.1 Hybrid Functionals ................................ ................................ ................... 99 6.3.2 Becke 3 parameter Lee Yang Parr (B3LYP) Hybrid Functional .............. 99 dependent (BMK) Hybrid Functional ................... 100 o B97 (B98) Hybrid Functional ....................... 102 6.4 Calculation Details ................................ ................................ .......................... 103 6.4.1 Gaussian 2 (G2) Method ................................ ................................ ....... 105 6.4.2 Gaussian 3 (G3) Method ................................ ................................ ....... 106 6.4.3 Complete Basis Set (CBS QB3) Method ................................ ............... 106
9 7 SURFACE MODIFICATION REACTION RESULTS ................................ ............. 109 7.1 Enthalpies of Formation ................................ ................................ .................. 109 7.2 Reaction Pathways ................................ ................................ ......................... 110 7.3 Summary ................................ ................................ ................................ ........ 112 8 CONCLUSIONS ................................ ................................ ................................ ... 117 LIST OF REFERENCES ................................ ................................ ............................. 120 BIOGRAPHICAL SKETCH ................................ ................................ .......................... 129
10 LIST OF TABLES Table page 7 1 Comparison of the enthalpy of formation at 298K for reactant and product structures ................................ ................................ ................................ .......... 113
11 LIST OF FIGURES Figure page 2 1 The structural domains of a kinesin motor ................................ .......................... 22 2 2 The structure of a microtubule filament ................................ ............................. 23 2 3 Components of a 13 protofilament microtubule filament ................................ ... 24 2 4 Hybrid design model approach for kinesin motors during transport .................... 25 2 5 Fluorescence microscopy images of gliding microtubules filaments .................. 26 2 6 Microtubule interactions within the Hybrid design model ................................ .... 27 3 1 Simulation tool flow chart for nanowire and nanospool formations ..................... 44 3 2 Fluorescence microscopy images of biotinylated microtubules .......................... 45 3 3 Schematic of microtubule movement on a hexagonal lattice .............................. 46 3 4 Initial, intermediate and final stage images in a simulation trial .......................... 47 3 5 The formation process of a nanowire ................................ ................................ 48 3 6 Angular and linear head to head microtubule intersections ................................ 49 3 7 The formation process of a nanospool ................................ ............................... 50 3 8 Confirmation of equal probability in lattice site visitations during movement ...... 51 4 1 Image of a simulated nanospool ................................ ................................ ......... 71 4 2 Images of initial, intermediate and final stages in a simulation trial .................... 72 4 3 Formation process of a simulated nanospool ................................ ..................... 73 4 4 Confirmation of microtubule movement with the turn probability values ............. 74 4 5 Analysis of output for 200 simulation trials on a 200 m x 200 m hexagonal grid surface ................................ ................................ ................................ ......... 75 4 6 Effect on the simulation output as the MT density is varied ................................ 76 4 7 Effect on the simulation output as the MT length is varied ................................ 77 4 8 Effect on the simulation output as the MT turn probability is varied .................... 78
12 4 9 Analysis of the number of nanospools generated versus time ........................... 79 4 10 Description of previous studies describing interactions proposed to initiate nanospool structures ................................ ................................ .......................... 80 4 11 Experimental images of nanospool formations ................................ ................... 81 4 12 Experimental size distribution of spool circumferences ................................ ...... 82 4 13 Comparison amongst experimental and simulation nanospool circumference length distributions ................................ ................................ .............................. 83 6 1 Flow chart of the iterative process to find a solution in the Hartree Fock Approximation ................................ ................................ ................................ ... 107 6 2 Flow chart of the iterative process to find a solution to the Kohn Sham Equations ................................ ................................ ................................ ......... 108 7 1 Optimized reactant structures ................................ ................................ ........... 114 7 2 Reaction pathway for thiophene and C 2 H forming a thiophene radical and C 2 H 2 ................................ ................................ ................................ .................. 115 7 3 Reaction pathway for thiophene and CH 2 forming a thiophene radical and CH 3 ................................ ................................ ................................ ................... 116
13 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONAL INVESTIGATIONS OF ORGANIC MATERIALS FOR HYBRID NANODEVICE AND OPTOELECTRONIC APPLICATIONS By Jasmine Davenport Crenshaw May 2011 Chair: Simon Phillpot Major: Materials Science and Engineering This dissertation examines two orga nic material systems, biotinylated microtubule filaments and thiophene. Biotinylated microtubule filaments partially coated with streptavidin and gliding on surface structures. We present a cellular automaton simulation tool that models the dynamics of microtubule gliding and interactions. In this method, each microtubule is composed of a head, body, and tail segments. The microtubule surface density, lengths, persist ence length, and modes of interaction are dictated by the user. The microtubules are randomly arranged and move across a hexagonal lattice surface with the direction of motion of the head segment being determined probabilistically: the body and tail segmen ts follow the path of the head. The analysis of the motion and interactions allow statistically meaningful data to be obtained regarding the number of generated spools, radial distribution in the distance between spools, and the average spool circumference lengths which can be compared to experimental results. This technique will aid in predictions of the formation process of nanowires and nanospools. Information
14 regarding the kinetics and microstructure of any system can be extracted through this tool by t he manipulation in the time and space dimensions Chemical reactions of thiophene with organic molecules are of interest to chemically modify thermally deposited coatings or thin films of conductive polymers. E nergy barriers are identified for reactive sys tems involving thiophene and small hydrocarbon radicals. The transition states for these reactive systems occurred through hydrogen abstraction. The results provide quantum mechanical level insights into the chemical processes that occur in the chemical mo dification processes described above, such as Surface Polymerization by Ion Assisted Deposition (SPIAD), electropolymerization, and ion beam deposition. Enthalpies of formation are calculated for organic molecules using B3LYP, BMK, and B98 hybrid functiona ls. G3 and CBS QB3 are used as standards in conjunction, due to their accurate thermochemistry parameters, with experimental values. The BMK functional proves to perform best with the selected organic molecules.
15 CHAPTER 1 GENERAL INTRODUCTIO N There is an enormous demand for devices to have enhanced properties at the micro and nanometer length scales. However, the smaller dimensions create difficult material challenges associated with high performance and very small scales. It is critical to understand the fundamentals of the structure of a material and its interactions in order to engineer properties for device usage. Thus, the field of Nan otechnology has emerged over the past few decades to produce new materials, devices, and instrumentation on the nanometer scale and develop ways to eliminate lim itations due to size reduction. The work described in this dissertation examines the fundamenta l properties of two organic materials systems. On the micrometer scale, the dynamics of microtubule filaments are elucidated. On the nanometer scale, the kinetics of chemical processes associated with the formation of thiophene in polythiophene thin films are examined. The human eukaryotic cell is composed of an intricate organization involving mechanisms for transportation to perform basic functions. Cytoskeletal filaments, such as microtubules, provide structural support to the cell and serve as pathways for the kinesin biomolecular motor protein to transport cargo such as organelles and vesicles during intracellular transport  Some of the fundamental processes associated with this intracellular transport system have been integrated into engineere d hybrid nanodevices. The key feature is the nanoscale functionality. Hybrid devices have been developed to drive self assembly; these incorporate proteins such as kinesin and microtubule filaments because they facilitate self assembly processes in biologi cal systems. In particular, by suitable device engineering, described in more detail in
16 Chapter 2, cargo such as mitochondria  endosomes  synaptic vesicles  and protein carriers  can be captured and transported to a destination via active transport and assemble nanostructures from a bottom up approach  Chemical reactions of thiophene wi th organic molecules occur during the production and chemical modification of thermally deposited coatings or thin films of conductive oligomers and polymers. One such material, polythiophene, is attracting much interest due to its stability at high temper atures, and because its molecular structure can be readily transfigured to achieve desired electronic properties  However, it is well known that thermally deposited polymeric thin films degrade in response to mechanical deformation and extreme thermal fluctuations [8 10] Chemical modification is one way to stabilize these films to external stresses without substantially altering their electronic and molecular structure. For example, modifications have been made to thiophene through surface polymerization by ion assisted deposition (SPIAD) by Hanley and co workers [11 15] The dissertation will use computer simulation methods to examine key aspects of these two systems. C omputer simulation s aim to mimic a few to several behaviors of a real system. Generated r esults can often compared to experimental results to obtain validation; moreover, the simulation can often explore phenomena or regions of parameter space that are not accessible to experiment, thereby extending understanding and providing predictions. Mathematical algorithms and key materials properties are embedded into the simulation to establish initial conditions. Computer simulations offer the capability of performing multiple tests on different length and time scales that may not be experiment ally accessible, or may be prohi bitively expensive to
17 probe. The amount of real time to complete a simulation vary in real time can range from a few seconds to days. Chapters 2 through 4 are dedicated to the study of microtubule filaments and their interactions to form wire and spool nan ostructures. Chapters 5 through 7 focus on surface interactions involving thiophene with organic radical C 2 H and CH 2 for the development of polythiophene thin films. Chapter 8 will integrate the findings in both project and project future work.
18 CHA PTER 2 KINESIN POWERED MOLECULAR SH UTTLE 2 .1 Introduction Motor proteins are enzymes that convert chemical energy into mechanical work. Chemical energy is received from the hydrolysis of adenosine triphosphate, ATP, and mechanical work is then generated. The three largest molecular motors are myosin, kinesin, and dynein. Myosin and kinesin are linear molecular motors that express mechanical energy as physical translation. Large concentrations of myosin motor proteins are found within muscle tissue of animals to perform contractions  Kinesin and dynein motors are both responsible for intracellular transport within the cytoplasm  Dynein motors are found in flagella used for cell motility, while kinesin motor s serves as actuators  to transport items such as organelles, RNA, and viruses. Kinesin motors are also involved in mitotic spindle formations and chromosome separations during cell division [1,19] The focus of this thesis will be kinesin motor proteins. Kinesin motors offer a number of unique capabilities to the field of Materials Science and Engineering. First, they can move unidirectionally along engineered tracks with their speed regulated by the concentration of adenosine triphosphate (ATP) fuel. Second, specific cargo can be captured from solution, attach to the motor protein, and thence transported  they can travel large distances over cytoskeleta l filaments called micr otubules without detachment. Structurally, conventional kinesin is composed of seven main components: two globular heads, dimerization parallel coiled coil dimer, flexible link domain, coil 1, kink
19 domain, coil 2, and tail domain (Fi gure 2 1)  T he movement of kinesin is driven by the globular heads. The two heads are structurally equivalent and walk in a hand over hand motion along the microtubule filament in 8 nm steps, the interval distance between tubulin dimers. The dimerization parallel coil ed coil dimer coils 1 and 2, and the kink domains are all coil regions with structural motifs of intertwined alpha helices. The link domain, also referred to as the neck, provides kinesin with flexibility during motion. The kink domain separates coil 1 fo rm coil 2. The cargo attaches to the motor tail during transport  Kinesin movement is driven by the hydrolysis of ATP. During hydrolysis, water is added to ATP causing the phosphate phosphate bonds to break yielding adenosine diphosphate (ADP), an inorganic phosphate (Pi), and 25 kT (100 x 10 21 J) of free energy. For each ATP molecule hydrolyzed, kinesin takes one step in movement. The globular head experiences a conformational change, and releases from the microtubule to attach to the next binding site on the microt ubule  2.2 T he Microtubule Protein Filament Microtubules are used as roadways for kinesin and dynein motors during transport and provide structural stability and support in eukaryotic cells. Microtubules have a pipe like configuration (Figure 2 2) with an inner diamet er of approximately 18 nm, an outer diameter of approximately 25 nm  and a length that can vary from a few nanometers to several microns [22,23] tubulin heterodimer subunits that bind together in a head to tail arrangement forming a protofilament [1,22] The protofi laments run parallel to the axis of the microtubule. A cylindrical formation results from the integration of protofilaments creating a sheet that wraps and binds to itself. The number of protofilaments per microtubule can vary from 8
20 to 19, but 13 is the m ost common number in cellular microtubules, as displayed in Figure 2 3  The tubulin dimer subunits within each microtubule are arranged in a slightly twisted hexagonal lattice, resulting in differing neighbor relationships among each subunit and its six nearest neighbors  Microtubules have two different ends positive end than at the neg ative end; as a result kinesin motors walk in the direction of the positive end  Taxol has been frequently used to stabilize microtubule lengths  2.3 Fluorescent Microscopy Fluorescent Microcopy (FM) is a visualization method that provides optical clarity and reve als the assembly, dynamics, and movement of organic and inorganic substances such as proteins. FM offers advantages in its use to specifically label individual or multiple molecules through the use of molecules such as fluorophore, green fluorescent protei n (GFP), and fluorescent beads. Photobleaching, a dynamic process where fluorescent molecules undergo photo induced chemical destruction upon exposure to excitation light and lose their ability to fluo resce, can plague FM  2.4 Hybrid Design Model The hybrid design approach (Figure 2 4) is a method for microtubule transport, which involves synthetic environments and kinesin motors. Kinesin motor proteins, which are adsorbed to an engineered surface, propel microtubules forward along fabricated tracks with a velocity of several hundred nanometers per second. Cargo, such as vesicles and organelles, can be attached to the microtubules at various locations using biotin linkers [2 5,29 37]
21 Functionalization of the microtubules with specific biotin linkers enables cargo loading  ; control of the ATP concentration governs shu ttle velocity  ; and tracks patterned on the surface select transportation paths  Functionalization of microtubules with biotin and partial coating with streptavidin enables cross linking of gliding microtubules after c ollisions (Figure 2 5); the surprising result of the self ultimately spools (Figure 2 6)  This result has broad significance since it shows that active transport by molecular motors can enable the acceleration of the self assembly process, the formation of non eq uilibrium structures, and the emergence of structure as a result of active transport processes 
22 Figure 2 1 The seven structural domains of a k inesin motor : two globular heads, the dimerization parallel coiled coil dimer, link flexible domain, coil 1, kink domain, coil 2, and a tail domain. Reproduced from Ref erence 
23 Figure 2 2. Structure of a microtubule and its subunits. (a) The and tubulin monomers come together to form a hetero dimer, which is the subunit of a protofilament. (b) Protofilament of a microtubule molecule. The plus end of the microtubule represents the direction of polymerization. (c) A cylindrical structure of microtu bule with approximately 13 protofilaments. Reproduced from Ref 
24 Figure 2 3. A 13 protofilament microtubule is positive and negative ends represent the direction in which tubulin subunits polymerize to the microtubule filament. The positive end polymerizes tubulin subunits at a faster rate than the negative end. Reproduced from Ref 
25 Figure 2 4. Hybrid model approach in which kinesin motors are attached to a surface and microtubules with cargo attached by biotin linkers are transported as the kinesin heads walk along the microtubules. Reproduced from Ref 
26 Figure 2 5. The assembly of gliding microtubules experimentally observed with fluorescence microscopy. At 10 seconds, the tip of one microtubule collides with the center of another microtubule and sticks. The tip is then forced to reorient and follow the other microt ubule at 20 seconds. This leads to the complete integration of the colliding microtubule into the leading microtubule or microtubule structure at 30 seconds. Reproduced from Ref erence 
27 Figure 2 6. Microtubule interactions within the Hybrid design model (a) Kinesin powered molecular shuttles rely on surface adsorbed kinesin motors to propel ca rgo carrying microtubules. Microtubules can be functionalized with biotin and streptavidin to enable cross linking between microtubules. (b) Collisions aggregates. Reproduced f rom Ref eference 
28 CHAPTER 3 CELLULAR AUTOMATON S IMULATION METHOD 3.1 Biomolecular Shuttle Transport Microtubules are observed to produce nanowires and nanospools from ATP driven kinesin motors in a non equilibrium state  Experimentally, transport and assembly processes of molecular shuttle motion face limitations. The positions of shuttles are identifi ed by fluorescence microscopy, wherein exposure of the fluorescently labeled microtubules to excitation light leads to photobleaching. A few dozen images, photobleaching  interferes with the smooth transport due to pho to induced cross linking of microtubules to the surface bound kinesin motors  A simulation tool that can encompass the spatial and time positions of microtubule shuttles in the absence photobleaching would be an asset. Details on the transport and self assembly format ions of microtubules into n anostructures can be extracted. This project involves the generation of a simulation tool to provide a fundamental understanding of the microtubule dynamics entailed in the self assembly formation of nanowires and nanospools. Thi s simulation tool will allow the following questions to be examined in order to extract information regarding microtubule dynamics: Can self assembly through the formation of binding interaction be simulated through computation to understand the fundamenta ls of nanowires and nanospools? Can computation produce comparable results with an experimental system? Can a depletion zone be present on the surface during nanospool formations? The simulation tool should allow us to capture different aspects of the expe rimental setup, such as a match in the surface dimensions to the microscope field of view dimensions, the number of microtubules (MT Density), the microtubule length
29 (average length being 5 m, and the microtubule persistence length ranging from 10 100 m [40,46,47] 3.2 Previous Computational Models Examining Microtubule Dynamics Shuttle trajectories on a planar surface have been experimentally characterized as worm like chain s with projection persistence lengths on the order of 10 100 m [40,46,47] To date, simulations have focused on the study of interactions between the gliding filament (microtubules or actin filaments) and the motors on the surface [48,49] ; the interaction of the gliding filament with a track defined by surface patterning  ; the interaction of the gliding filament with external forces  3.2.1 Fluorescent Microscopy Fluorescent Microcopy (FM) is a visualizat ion method that provides optical clarity and reveals the assembly, dynamics, and movement of organic and inorganic substances such as proteins. FM offers advantages in its use to specifically label individual or multiple molecules through the use of molecu les such as fluorophore (green fluorescent protein (GFP) and fluorescent beads. Photobleaching, a dynamic process where fluorescent molecules undergo photo induced chemical destruction upon exposure to excitation light and lose their ability to fluoresce, can plague FM  3.2.2 Off lattice Monte Carlo Simulation Nitta et al uses an Off lattice Monte Carlo simulation to model the gliding motion of molecular shuttles involving microtubules and kinesin motors  In the 2 D simulation code i s written in Fortran 90. The forward direction of motion the molecular shuttle is determined by random fluctuations at each time step. Experimental parameters are incorporated into the simulation in order to embed information regarding
30 the microtubule and kinesin structures. The parameters included are the time averaged velocity, v avg persistence length, L p motional diffusion coefficient, D v flu The simulation parameter values of v avg D v flu and L p respectively, 0.85 m s 1 2.0 x 10 3 m 2 s 1 and 11 1 m  The time step in the off lattice Monte Carlo simulation m odel is set to 100 ms, which is the time a kinesin motor needs to take 10 elementary steps of 80 nm. The movement of a microtubule is determined from the trajectory of its leading tip. The microtubule tip trajectory has a fixed persistence length and is ba sed on a random walk in the Monte Carlo Simulation. During each time step, the microtubule step distance, r, random fluctuations. The mean and variances of the step distances and angular changes are displayed in E quations 3 1 through 3 4  = v avg (3 1) = 2D v flu (3 2) = 0 (3 3) (3 4) This model reproduces properties of experimentally realized systems. The molecular shuttle movement is visualized as whole networks, where a rotational design of individual guiding structures is generated. However, this simulation model neglects interactions between microtubules.
31 3.3 Cellular Automaton 3.3.1 Background on C ellular Automaton The Cellular Aut omaton (CA) paradigm is a simulation method that examines spatial interactions amongst structures moving over evenly spaced, discrete cells. The lattice type can be varied based on the system of interest. A few examples of the lattice types include: square triangular, and hexagonal. A cellular automaton requires a lattice of adjacent cells covering a portion of a one or multi dimensional space; a set of variables attached to each lattice site and the local state of each cell at each discrete time step; the state of each cell, location, and cell neighbors are all governed by user supplied transition rules to observe local interactions  The first universal cellular automaton was offered in the work of von Neumann, where a 2 dimensional universal and self reproducing CA simulation was used with 29 states to examine self reproduction  The CA approach provides advantages, including low computational cost; simple implementation  ; the ability to simulate large system sizes  ; and the ability to follow the evolution in time over experimentally accessible periods  The cost paid for this efficiency and simplicity is that CAs incorporates the underlying physical phenomena through simple rules. CAs provide insight into the fundamental physical features o f the system from the microscopic and macroscopic level. There are limitations within the Cellular Automaton method. Whether the rules transition rules defined do or do not capture the essential physics has to be determined a posteriori by the degree to wh ich the simulation results obtained match experiment or results from more sophisticated simulation methodologies  In CA, the models are purely kinematics and no flow dynamics are involved 
32 3.3.2 Types o f CA methods There are a number of Cellular Automaton methods that have been developed to study various systems and perform operations. For example, the lattice gas cellular automata (LGA), Frisch, Hasslacher and Pomeau (FHP) automata, and Lattice Boltzman n Automata (LBA) are all used to examine fluid flow and dynamics. The Lattice Gas Cellular Automata (LGA or LGCA) method simulates fluid flow by solving the Navier Stokes E quations  In this model, the lattice sites are associated with different particle states. Each state is characterized by particle velocities. The state of each site is determined with logic with consideration taken amongst each site and its neighbors. Propagation and collision are then ca lculated during each time step. The Frisch, Hasslacher and Pomeau (FHP) automata was introduced is to examine binary particles propagating through a hexagonal Bravais lattice at discrete time ste ps and a given unit velocity  During propagation, all particle movement is based on the direction of its velocity vector and the lattice connection lines between the six nearest neighbors. Particle mass and momentum are conserved during collisions. In FHP, an exclusion principle is im posed that does not allow the location of more than one particle per cell at the same time. Macroscopic quantities, such as pressure and velocity, can be extracted from particle distributions, appropriate time, and space averaging procedures  In Lattice Boltzman Automata (LBA) is similar to FHP automata, where fluid dynamics are examined. However, the two methods differentiate in that the LBA method observes fluids over long period time and large areas. In LBA, computations are performed on smaller grids with less iteration. The density distribution in the fluid flow
33 does not chan ge, and the final density can be obtained without any time and space averaging processes  3.4 F ortran 90 The programming language used to write the simulation is FORTRAN 90. FORTRAN 90 offers advantages such as: array features to keep track of spatial and time features in Cellular Automaton, modules, and flexibility in its output to be transferred into other programming and visualization environments. Fortran 90 contains array through the use of ALLOCATE and DEALLOCATE statements. This feature allows efficient storage allocation by issuing space according to the size of the problem at hand within the code  Pointers were also added to Fortran 90 to allow data objects to be declared with an attribute so that an object does not have any storage until an existing target object. In the case of an array, only the rank is declared initially and a shape is acquired when it is associated with a target  Modules are collections of data, type definitions, and procedure definitions. Modules introduce a course of action for adding intervals. An interface block communicates with the compiler to associate a functi on with an operator, and indicates similar code for doing other operations on interval data types, and similar code for operations between reals and interval data types. This addition allows flexibility for complicated mathematical equations to be incorpor ated into the code for calculations  large character sets with more than two precisions for real and complex values. Unlike Fortran 77, a
34 complex number must be supported with the same set of precisions as a real number  3.5 Development of the Simulation Tool 3.5.1 Overview The simulation tool developed here presented is influenced by the Off lattice Monte Carlo simulation from Nitta et al In particular, it incorporates interactions amongst microtubules in order to understand the self assembly process in the formation of nanostructures. A flow chart of the simulation components are displayed in Figure 3 1. 220.127.116.11 Definition of lattice A hexagonal lattice is used in our simulations as it provides a good compromise between simplicity and the ability to realistically describe dynamical movement of the microtubules, including their ability to twist and turn over the surface, and to form comp lex structures including spools. The user provides information that defines both the initial conditions for the simulation and the rules under which the CA evolves. In terms of initial conditions, the user inputs the dimensions of the hexagonal lattice, th e number of microtubule chains and the length of the microtubules. For the proof of principle simulations discussed in this chapter, the total dimensions of the system selected are 80 m x 80 m to match the field of view in the FM microscope and 200 m x 80 m The total simulation time is chosen to be similar to that accessible to experiment: typically of the order of minutes to a few hours. In the version of the code described here, all of the microtubules will either have the same length or a modificat ion to include a distribution in lengths. Details regarding the microtubules will be discussed further in Chapter 5.
35 The lattice unit size is also chosen to be similar to the experimental length scales. The center to center distance between two interacting biotin/streptavidin functionalized microtubules is approximately 35 nm  The distance between parallel lines of the lattice could therefore reasonably chosen to be 70 nm, so that two parallel microtubules with a center to center distance above 35 nm should be represented on different lattice points. The lattice unit, lu, size of the hexagonal lattice is actually chosen as 80 nm, which corresponds to the pixel dimensions of the experimental images with the 100x objective most usually used (Fig ure 3 2). Periodic boundary conditions are applied to the hexagonal lattice surface to generate a continuous surface of movement for the microtubules. The lattice is hexagonal, but the simulation surface cell is actuall y a rhombus (Figure 3 3). As a microtubule reaches the edge of the surface it re enters on the opposite side moving in the same direction. 18.104.22.168 Simulation parameters The simulation tool incorporates parameters that mimic the experimental setup to observ e cellular transport and microtubule dynamics. The simulation parameters that can be controlled are the dimensions of a hexagonal surface, the number of microtubules, the number of elements in an individual microtubule, the turn probability (which correspo nds to the microtubule persistence length), and the number of time steps. Experimentally, forces from the kinesin motors constrain the microtubule assemblies to remain metastable in the presence of ATP and stable in its absence  The assembled protein structures, nanowires and nanospools, are stabilized to maintain their mesoscopic conformation even after drying due to covalent cross linking  The driving forces present within the system are from the kinesin motor and the short range
36 interactions of biotin streptavidin amongst the microtubules. The kinesin motors provide the motility of the microtubules on a surface, where the biotin streptavidin bonding drives the nanowire and nanospool forma tions. The effect on the driving force of the kinesin motor is represented by the motion velocity of microtubules in the simulation. It is found experimentally that microtubules are propelled in movement by the kinesin motor at a velocity ranging on an ave rage o f 800 nm/s [61 63] In the simulation, each microtubule is set to move throughout the hexagonal surface at a rate of 800 nm/s. The distance between each lattice site is 80nm and thus each time step corresponds to 0.1 second of real time. The tim e step is defined by the ratio of lattice size and microtubule velocity. The microtubule velocity of 0.8 m s 1 is selected from the microtubule velocity distributions of Nitta et al  and is set as a standard in our simulation tool. Therefore, each simulation time step corresponds to 0.1 seconds, during which all indep endent microtubules move by one lattice site. The system size of 200 m is an appropriate length scale for the surface dimensions, because it enables the formation of several spools per run, which reduces boundary effects. The effect of the streptavidin bi otin bonds is initiated when two or more microtubules intersect and occupy the same lattice position during the same time step. The bond formation allows independent microtubules to merge together at the point of intersection to form an elongated nanowire, composed of multiple microtubules. Given that the ratio between the microtubule trajectory persistence length of 0.1 mm  and the lattice size of 80 nm is more than 1000, the probability of a 180 o revers e in direction within one step can be neglected on the time scale of the experiment.
37 22.214.171.124 Description of simulated microtubules Experimentally, the average microtubule length is 5 m  which corresponds to a few dozen subunits on the 80 nm lattice. In the simulated microtubule, the first subunit is defined as the head subunit, the last subunit is the tail subunit; all subunits in between are labeled as the body subunits. Each microtubul e has an orientation on the lattice, distinguishing the head from the tail, as shown in Figure 3 3. The direction of motion of each microtubule is determined by the head subunit. To set up the initial arrangement on the hexagonal lattice a random generator selects one of six orientations for each microtubule, as shown in Figure 3 3. In the initial simulation set up, the microtubules are not allowed to overlap, cross or intersect. Initially, the microtubules are randomly arranged and are aligned straight bas ed on their orientation on the surface as in Figure 3 4. Since shuttle trajectories on a planar surface have been experimentally characterized as worm like chains with projection persistence lengths on the order of 10 100 m [40,46,47] the parameters necessary to create trajectory ensembles with realistic average behavior and fluctuations around it are available. 126.96.36.199 Motion of microtubules The most important inputs are the ru les under which the cellular automaton evolves. The rules implemented are designed to match the fundamental mechanisms inferred from experimental observations. If the simulations then reproduce the fundamental of the experiments it can be assumed, though n ot guaranteed, that the rules do indeed reflect some aspects of the underlying physical processes. If the simulation results depart significantly from the experimental results, then it can be
38 inferred that either one of more of the rules does not reflect t he underlying physics, or that a physically important phenomenon has not been included in the rules, or both. Periodic boundary conditions are applied to the hexagonal lattice to generate a continuous surface of movement for the microtubules. As a microtub ule reaches the edge of the hexagonal lattice it re enters on the opposite side moving in the same direction. The user has to define two additional conditions in order for movement to begin: the turn probability for each microtubule and the time duration o f the simulation, i.e. the number of steps over which the CA will evolve. 188.8.131.52 Turning of microtubules Microtubules have random movement over a surface, but the movement is also governed by its persistence length in experimental observations. In the sim ulation, the turn probability accounts for the random motion of microtubules through the number of turns made by the microtubule by incorporating the microtubule persistence length. Limitations are present with the turn probability parameter because it is based off an approximation of the microtubule persistence length. The turn probability, p, is calculated from considering the cosine of the average deviation from the original direction of after 1 step by the microtubule filament. Given the relationship between the experimental persistence length to the turn probability in E quation 3 5 the definition of the persistence length is displayed in E quation 3 6 where is the lattice unit length. cos< > = 1 2p+2pcos(60) = 1 p (3 5) exp( l/2L p ) = (3 6) p = /2L p (3 7)
39 As a result, the probability of the microtubule making a sixty degree turn can be calculated in terms of the persistence length and latt ice unit parameter in E quation 3 7. Therefore, for the simulations provided the turn probability is set to 0.0005 meaning in the microtubule changes to make a 60 degree turn in direction approximately once every 2,000 steps, and a 120 degree turn once ev ery 4,000,000 steps. 3.5.2 Binding Interactions Transition rules are implemented into t he simulation tool to account for interactions of between microtubules that are dependent and independent from each other. A microtubule is labeled as dependent if its tip intersects with the body subunit of a second microtubule. The microtubule that is in tersected is labeled as independent. If subunits of two independent microtubules occupy the same lattice position at the end of a time step, then a flag is generated indicating that a chain intersection has occurred. In the current implementation, based on experimental observations, it is assumed that there is a 100% sticking probability on such an intersection; this could easily be changed in the models. On intersection, the chain whose head intersects the body of fashion as shown in Figure 3 5 after which they act as a single elongated chain, with the leader chain determining the direction of movement of both chains. As the microtubules merge, there is an overlap, where the leader and follower microtubules occupy the same lattice positions. This process can occur multiple times such that a single leader chain may have a large number of parti ally overlapping follower chains.
40 If the head subunits of two independent microtubules meet at the same site, they cannot both become the leader chain. Therefore, a special rule is invoked in the case of head to head collisions of leading chains. If the he ad to head collisions occur at an angle (60 o 120 o 240 o or 300 o ), an additional rule is applied to handle the merger of two microtubules into a single elongated chain. In this case, one of the two microtubules is randomly selected to become the leader chain as illustrated in Figure 3 6 A by the purple chain. The follower chain, in cyan, is dependent on the leader chain to define its movement. If a linear head to head collision is confirmed, a random generator selects a new direction for one intersecting microtubules chosen at random, in order for it to prevent the collision (Figure 3 6 B ). 3.5.3 Treatment of Nanospool F ormation s The experimental results of Hess et al  have shown that nanospools can form when the tip of a microtubule chain intersects with a part of the same chain. In order to mimic the same behavior in the CA simulation, we have cre ated a corresponding rule: nanospool formation is initiated when the head subunit of a microtubule intersects the body or tail subunit of its own chain, or a subunit of one of its attached follower chains. In Figure 3 7, the head subunit corresponds to the point of intersection where the microtubule merges to begin to form a nanospool. The subsequent evolution involves the tail and follower chains wrapping up to form the spool; the total perimeter of the spool is then determined to give a measure of the siz e of the microtubule chain. The simulation ends when either the prescribed number of steps is reached, or when all of the chains have formed into nanospools. At the end of the simulation, we determine the number of nanospools that have formed and determine the circumference of each spool. Multiple simulations are performed under the same rule set, but with
41 different random initial conditions in order to obtain good statistics. The key quantities tracked in each simulation and analyzed statistically are the number of spools generated, the spool length, and the radial distribution of the spool locations. 3.6 Simulation Memory Consumption As mentioned above, one of the key appeals of the CA method is its computational efficiency. In particular, as demonstrated here it can access both the length and time scales of the experiments. The greatest memory usage in the simulation is for the two dimensional arrays, which are used to track the lattice positions of each microtubule and the dimensions of the lattice surfa ce. The memory usage scales linearly with the microtubule density, microtubule length, and linear dimension of the lattice. Improvements can be made in the performance of the next generation code by conserving memory. Currently, 2 D arrays are defined to p erform specific tasks and calculations in the simulation. Ultimately, memory space is consumed while the arrays are no longer in use for the remainder of the simulation. Local, temporary arrays can be implemented to reduce memory consumption. The CPU time was found to scale quadratic with the number of microtubules, microtubule length, and linear dimensions of the hexagonal lattice. 3.7 CA Algorithm Verification There are a number of simple indications to verify the correct operation of the simulation tool. The indications include: the simulation will begin with a set number of microtubules and end with nanospools; microtubules are unable to take on backwards motion onto itself; independent microtubules cannot experience head on head interactions during move ment; the number of intersections should not be more than the number of microtubules.
42 The simulation begins with a set a number of microtubules randomly arranged over a hexagonal lattice. As time evolves, the microtubules will randomly move throughout the lattice to encounter binding interactions with other microtubules. By defining long simulation time periods, an appropriate time will be issued for all the microtubules interact and form nanowires and then nanospools. During movement, microtubules are unab le to move in the backwards direction in movement. Backward motion consequently forms interactions that are not of interest in the scope of this project. Therefore, if 180 o degree turns occur in the simulation, I know that the motion indicated are not of i nterest in this project. Decision points are embedded in the simulation to generate flags that indicate head to head collisions. If a head to head collision occurs during movement, a clear indication is generated that a problem exists in the current trans ition rules for movement. The head to head collision represents a movement violation in to identify all microtubules involved during binding interactions. If each microtubule is unable to be properly identified, then the current rules employed cannot accu rately account for all spatial interactions of microtubules during self assembly. A count in the number of simulation interactions can provide insight in verification of the simulation. If the total number of interactions is greater than the initial number of microtubules, a problem exists in the evolution of the simulation. Within the simulation, each microtubule is to independently interact with the one another. If the number of interactions is incorrect, a flag is generated in the code and the simulation stops. The implementation of this rule helps to account for the number of microtubules in space and time.
43 There are key measures within the simulation to confirm if the simulation is operating in the correct manner. During the simulation, each lattice pos ition should have the opportunity of being occupied b y a microtubule filament. In Figure 3 8, one microtubule is randomly placed on an 8 m x 8 m (100 lu x 100 lu) hexagonal surface. The position of the head subunit is tracked for its occupancy on the lat tice surface. A counter is placed at each lattice position to tally the number of visitations. Figure 3 8 verifies equal probability for lattice coverage at an average of 1100 occurrences over the entire lattice positions.
44 Figure 3 1. Flow chart of the simulation to explain the spatial calculations that account for nanowire and nanospool formations amongst the microtubules.
45 Figure 3 2. Fluorescence microscopy images of biotinylated microtubules, partially coated with streptavidin and gliding on surface adhered kinesin motor proteins form linear nanowires and circular nanospool structures. Reproduced from Ref erence 
46 Figure 3 3. Schematic of microtubule movement on the hexagonal lattice. Each simulated microtubule is composed of a head, a tail, and number of body subunits. The calculation of new lattice coordinates during movement is based on the orientation of the microtubule. A weighted random number determines the movement for the micr otubule; after the head subunit moves each subsequent body and tail subunits will move into the position of the pre ceding subunit.
47 Figure 3 4. Initial, intermediate and final stage images in a simulation trial. Initially, two hundred microtubule filaments are randomly arranged on an 80 m x 80 m hexagonal lattice surface. Nanowire and nanospool formations are detected in time step 350 (35 seconds). Lastly, the simulation comes to a completion with nanospool str uctures remaining in step 50,000 (5,000 seconds).
48 Figure 3 5. The formation process of a nanowire. At time step 15 (1.5 seconds) two independent microtubules (yellow and purple) intersect. The head subunit of each microtubule occupies the same lattice position. The purple microtubule becomes the follower, and the yellow microtubule becomes the leader. At time step 15, the tip of the purple microtubule merges into the yellow microtubule from the point of intersection. At time step 45 (4.5 seconds) the positions of the two micr otubules overlap and occupy the same lattice positions (the overlapping subunits form the purple microtubule is not shown for visibility reasons). In time step 75 (7.5 seconds) the two microtubules have fully merged into an elongated nanowire.
49 Figure 3 6. The angular and linear head to head microtubule intersections are special cases in the assembly process. (a) In angular intersections two independent microtubules approach in an angular orientation in time step 6 (0.6 seconds); in time step 8 (0.8 sec onds), intersection occurs and is detected; one of the two microtubules is randomly determined to lead in the movement as seen in time step 12 (1.2 seconds); the two microtubules then merge in time step 39 (3.9 seconds). (b) In linear head to head intersec tions two independent microtubules approach in a linear orientation in time step 6 (0.6 seconds) and intersect in time step 8 (0.8 seconds). The intersecting microtubule will select a new direction for movement in order to avoid intersection. Normal inters ection will end up resulting because the second microtubule will intersect the microtubule that has changed course in movement in time step 35 (3.5 seconds).
50 Figure 3 7. The formation process of a nanospool: At time step 10 (1 second) an independent n anowire is approaching to cross itself. At time step 12 (1.2 seconds) the nanowire intersects itself to initiate the nanospool formation. The lattice position of the nanowire head subunit of the nanowire will now take on the lattice position of the body su bunit it intersected. As time evolves, the head subunit will overlap itself to form a nanospool at time step 47 (4.7 seconds)
51 Figure 3 8. Confirmation of equal probability in lattice site visitations during movement A microtubule is randomly placed on a 8 m x m (100 lu x 100 lu) hexagonal lattice surface. A counter is used track the positions of the microtubule head subunit and determine the number of lattice visitiations on the surface during movement. Over a span of 1,000,000 seconds (10,000,000 time steps), the microtubule equally visits all 10,000 lu positions at an average of 1100 occurrences.
52 CHAPTER 4 MT SIMULATION RESULT S AND DISSCUSSION In C hapter 3 a simulation tool has been developed to capture the mobility of microtubules and the binding effects of microtubules during the self assembly of nanowires and during nanospool formations. These simulations are carried out on a two dimensional hexagonal l attice. The initial location of each microtubule is defined by a position and a head to tail orientation determined at random; each microtubule then moves over the lattice, encounters other microtubules, and through interactions described by specific rules finally b ecomes part of a nanospool (Figure 3 7 ). In this c hapter, three scenarios are analyzed. The first scenario involves simple systems and situations designed to validate the simulation methodology (Section 4.1). The second scenario increases in com plexity by examining the effects of system variables, including the microtubule density, length, and trajectory persistence length (turn probability) on the self assembly of nanostructure formations (Section 4.2). Finally, the results of the simulations ar e confronted with experimental data in the third scenario (Sec. 4.3). Lessons will be learned regarding the strengths and weaknesses of the simulation tool, and the information that can be extracted from the simulation tool. Predictions concerning the self assembly will also be identified by systematically changing the system variables. Lastly, comparisons will be made amongst simulation and experimental results to test the capability of the tool. Microtubules with biotin linkers, and partially coated w ith streptavidin, can crosslink after collisions; the surprising result of the self shuttles is the formation of extended wires and ultimately spools (Fig ure 3 2 ). The cause of the initiation of the nanospool formatio ns is unknown; however, experimental
53 observations have led to a number of theoretical suggestions. Microtubules have been 750 m along a microtubule filament  A microtubule tip can get stuck on the surface at a nonfunctioning kinesin motor that is unable to transport during pinning [64,65] ; Brownian bending has also been detected, in which thermal fluctuations in the microtubule tip are measured to account for the micr otubule trajectory as a worm like chain  Simultaneous intersections have been recognized in spool formations, where multiple microtubules intersect and crosslink forming a triangular contour that relaxes into a ring like structure  4.1 Scenario 1: Validation of the S imulation T ool In the simulation, the microtubule movement and interactions generate non circular nanospools w ith complex configurations (Figures 4 1 and 4 2). The simulated nanospools produced do not take on ring shapes as in experimental FM images. The presence of the he xagonal lattice cannot support an entirely smooth ring. More importantly, in the simulation the microtubules have not intrinsic stiffness. As a result, there is no tendency for the irregular shaped closed loops to be smooth out into rings. The simulated na nospool formation process allows the random movement trajectory of a microtubule that intersects itself. The head subunit of the microtubule corresponds to the point of intersection where the microtubule merges to begin to form a nanospool. The subsequent evolution involves the tail and follower chains wrapping up to form the spool; the total perimeter of the spool is then determined to give a measure of the start of the microtubule chain (Figure 4 3).
54 Hundreds of microtubule chains are encompassed within a nanospool, produced from intersections with the leading microtubule during the formation of a nanowire or with the nanospool itself. In the simulation, all of the constituent microtubules occupy the same lattice coordinates; however, not all of the nanosp ools are reflected in the simulated nanospool images. In order to achieve an optimal visual appearance and clarity, only regions of overlap are displayed in Figure 4 1. The formation process from the simulation is consistent with the development process ob served in experiment. The main difference between the simulation and experimental systems is that the experimental nanospoo ls tend to be circular (see Figure 3 2). 4.1.1 Persistence Length/Turn Probability This section will confirm how the turn probability parameter in the simulations corresponds with the experimental setup of Hess et al Experimentally, challenges are presented when determining the exact probabilities of 60 and 120 turns of a microtubule within an 80nm step due to the few occurrences. In stead, small directional change distances from the microtubule tip are measured for 60 and 120 turns. Equation (4 1) displays how the directional changes are correlated with the persistence length. Consequently, the turn probability parameter, p60, for 6 0 turns can be ca lculated using Equations 4 2 and 4 3. ~ 1 (4 1) (4 2) ~ (4 3)
55 The probability for 120 turns is simply the square of the probability of 60 turns (Equation 4 4) : (4 4) From analysis of FM images yield a persistence length of 90 m. The 60 turn probability for the microtubules was calculated to be 0.00045 (0.045%), where the microt ubules were observed to be transported over a kinesin coated surface at a rate of 800 nm/s (v avg ) [61 63] The persistence length of the microtubules falls lower than the persistence length of 111 m  because of the presence of defective kinesin motors on the surface. Defective kinesin motors are able to bind b ut not transport a microtubule. 4.1.2 Valida tion of Nanospool Formation Changes in path trajectories also confirm the correct functioning of the simulation. The product of the number of microtubules and the number of time steps should equal the total number of trajectory changes (60 o and 120 o turns and the forward direction motion) on the lattice in E quation (4 5). (Lattice movement occurrences) = (# of microtubules)*(# of time steps) (4 5) The sum of the sixty or one hundred twenty degree turns divided by the total number of lattice movement occurrence should equate to the initial turn probability values set by the user. Figure 4 4 shows an example confirming the microtubule movement with the 60 o and 120 o turn probability parameters set in the initial conditions. 4.1. 3 Number of Spools Generat ed and Nanospool Circumferences Upon completion of a simulation run, the surface is littered with nanospools. In order to gain a better understanding of the self assembly process, the number of spools
56 formed is evaluated. Meaningful statistics are obtained through the repetition of the simulation for 200 trials. Figure 4 5 A shows the resulting histogram for the total number of spools produced in each trial. Each trial has an initial condition of 1250 microtubules of initial length of 5 m, a turn probabilit y of 0.045% (persistence length of 90 m), and runs for a simulation time of 10,00 0 seconds (100,000 time steps). The probability distribution in the number of spools produced on the simulated surface can be fit with a Poisson distribution. The Poisson dis tribution is selected for use due to the following conditions: spool formation is independent of kinesin; the probability to have a nanospool for each surface location is small; the number of spools produced during any simulation trial of 10,000 seconds is independent from every other trial; The average number of nanospools generated is used to characterize the Poisson distribution over a specified region. If the average number of nanospools per area is the average number of nanospools generated in an ar ea of size A is Therefore, the probability that at least one nanospool forms on the su rface is given in Equation 4 6. p A A ))] (4 6) The results in the number of spools were fit with a Poisson distribution to confirm its trend. Fi gure 4 5 A fits well with the Poisson distribution and the average number of nanospools generated corresponds with the peak of the curve to be 12. Figure 4 5 B displays a histogram of the output generated from simulated nanospool circumferences. While most s pools have a circumference of less than 100 m, a small peak is present at circumferences near 200 m. These large spools correspond to straight wires which have grown very long and eventually intersect with themselves through the periodic boundary conditi ons. Due to the clear separation in the
57 spool circumferences, spools generated from periodic boundary conditions can be readily identified and discarded from further analysis. The size of the generated spools are consistent with experimental results [41,66,67] where the average spool circumference is 10 m  The effects of change parameters, such as the number of microtubules, turn probability, and the microtubule length, will be discussed below in s ection 4.2 (s cenario 2 )  4.1.4 Nanospool C ircumference Radial Distribution The relative locations of spools are analyzed inside the simulated surface by c omputing the radial distribution function, g(r), in 2 D. The g(r) provides information on the microstructure of a material and verifies the presence or absence of order in the location of nanospool formations. For these simulations, the radial distribution function is defined as the ratio of the density of spools at separation r to the density if the spools were uniformly distributed. Thus, a value of g(r) = 1 for all values of r corresponds to a random distribution (as in an ideal gas); deviations from all g(r) = 1 signify either a local deficit in neighbors, g(r) < 1, or a local excess, g(r) > 1  The bar chart in Fig ure 4 5 C compared t around each spool, formed by all the microtubules in the vicinity being swept up into a single nanospool. The plateau region in the radial distribution function for the spools fo rmed, as determined at the end of the simulations (Figure 4 5 C ) indicates the absence of long range order among the nanospools; the average g(r) over the range r = 35 100 m is 0.99, indicating a random distribution and ideal gas like behavior.
58 4.2 Scenario 2: Effects of Simulation Parameters on Nanospool Formations The microtubule CA tool is used to characterize the spool formation process and to identify emergent properties of the self assembly process for microtubule nanowires and nanospools. The parameters of interest in the simulation include the initial microtubule length, microtubule density, and the microtubule persistence length (turn probability). In experiment, the microtubule density and length parameters can be controlled through the init ial concentrations of tubulin and taxol. The microtubule persistence length appears to be an intrinsic property of the system and is typically extracted from AFM images of the microtubules by analyzing the contour of the filament or the end to end distance ; however, the persistence length can be altered by electrostatic repulsion between charged monomers  The influence of these parameters will be studied for their effects on the spool formation process through the quantification of the number of spools, the circumference length of the spools, and the distribution of spo ol locations as a function of time. Upper limits of the efficiency of the self assembly process will be established in the simulation. These will be used to provide insights as to the presence or absence of experimental defects in the kinesin coverage of t he engineered surfaces; the dimensions of the largest assembled nanospool will be determined. The simulations will have 625, 1250, and 2500 microtubules randomly arranged on the hexagonal surface. The number of microtubules set in the simulation correspond s with the number of microtubules within the field of view of fluorescent microscopy images. Various trajectory persistence lengths will be explored; these will be characterized by 60 turn probabilities of 0.0%, 0.045%, and 4.5%. Within the simulations, m icrotubule filaments are set to have an initial length of 2.5 m, 5 m, and
59 10 m. These three lengths were selected out of a range in the distribution of microtubule lengths (0.8 m 23 m) observed in the experiments of Jeune Smith et al  The system size of 200 m is selected as an appropriate length scale for the surface dimensions, because it enables the formation of several spo ols per run, which reduces boundary effects. 4.2.1 Simulation Parameters The effects of the parameters on the nanospool format ions are reported in Figures 4 6 through 4 8. The number of generated nanospools is shown to correlate with microtubule interactio ns. As the microtubule density and length are increased, the microtubules are more prone for collisions (Fig ures 4 6 A and 4 7 A ). Increased turn probabilities allow more turns to occur on the lattice and ultimately self intersection or intersectio ns with ot her microtubules (Figure 4 8 A ). The MT length parameter greatly impacted the spool circumference, where longer microtubules are proportional to increased nanospool circumferences (Fig ure 4 7 B ). The size of the generated nanospools generated from microtubul es with an initial length of 5 m is discovered to consistent with experimental results [41,66,67] where the average spool circumference is 10 m  Figure 4 8 B confirms that decreased turn probability percentages yield longer nanospools. More specifically, ~95% of all the nanospools generated for turn probably 4.5% have a spool circ umference of 2 m. In observing simulation output movies, the increased turn probability allows for self intersections where the microtubules have an increased probability of making 60 o turns on the lattice yielding small circumference nanospools. In Figur e 4 6 B the microtubule
60 density is shown not to have an impact on the nanospool circumference due to the small variations in the spool circumference frequencies. The relative locations of spools can be analyzed by computing the radial distribution function g(r), in 2 D, which is defined as the ratio of the density of spools at separation r to the density if the spools were uniformly distributed. All of the radial distribution function curves in Figures 4 6 C 4 7 C and 4 8 C yield a random order at long radi al distances where g(r) = 1. However, trends are identified to characterize the depletion zones. Increased microtubule density, length, and turn probability parameters all generated smaller depletion zones. In Figure 4 7 C the radial distribution function curves for 5 m and 10 m nearly overlap. A difference of 10 m exists amongst their depletion zones, but after a radial distance of 40 m the radial distribution function curves overlap. 4.2.2 Nanospool Formation Versus Time The number of spool formations generated as a function of time is investigated in Figure s 4 9 A and 4 9B The simulation begins with no interactions (zero spools) and the number of spools quickly increases and level off to the total number of nanospools. The majority of the nanospools are formed within the first 60 seconds as a result of multiple interaction s on the lattice. However, there are cases where a few microtubules have not formed nanospools after one minute and addition time is needed for a microtubule to intersect a nanospool on the surface. The trend in the number of spool formations as a function of time is comparable to Langmuir adsorption isotherms, where the fraction of surface sites of adsorbed solute molecules are measured as a function of molar
61 concentration  In Figure 4 9 B a logistic function is determined to fit well with the spool formation trend. 4.3 Scenario 3: Application of the S imulation T ool in the Dependence on Trajectory Persistence Length The analysis below is a joint experimental and simulation program  The results are presented in an integrated manner so as to be most effective. The experiments discussed were performed in the group of Prof. H. Hess. 4.3.1 Introduction science  One approach to overcome the limitations in component size, assembly speed and structural characteristics of chemical self assembly is to utilize active transport rather than diffusion as the mechanism to achieve the recruitment and assembly of building blocks into a larger structure  Active transport, for example driven by molecular motors, can move larger components faster than diffusion. Also, the spectrum of forces exerted during and after assembly by active transport is dramatically different from the spectrum of thermal forces present during and after equilibrium self assembly. As a result, the assembly of non equilibrium structures, such as structures under high internal strain, should be possible. A simple model system for active self [41,73] Biotinylated microtubules form w hen thousands of biotinylated tubulin proteins polymerize into tubular filaments with a diameter of 25 nm and a length of several micrometers, and they can bind to kinesin motor proteins adhered to a surface 
62 because biotin streptavidin biotin cross links can form. When gliding sticky microtubules collide, they assemble into elongated bund les. When the tip of a bundle encounters the middle of the bundle after a sharp turn, 4 1 0 ). These spools typically have diameters of a few micrometers, which means that their formation from bundles of microtubules, each having a persistence length on the order of millimeters [74 76] consumes a significant amount of energy. This energy amount on the order of thousands of kT has to be supplied by the hydrolysis of ATP and transduced by the motors. Spool formation is thus a consuming self  since the spool is stable ev en when the energy flow ceases. The first question of interest here is if the initiation of spool formation is the result of thermal fluctuations or the result of motor action. The answer is of considerable interest, because compared to a process purely co ntrolled by the motor action an active self assembly process which utilizes a thermally activated process as a rate limiting step is significantly less amenable to engineering control. The second question we seek to answer is: what controls the size of the spool formed? To answer these questions, we measure the size distribution of microtubule nanospools, and show, using computer simulations of gliding bundles, that thermal fluctuations in gliding directions result in spools which are four times larger in d iameter than the experimentally observed spools. We then discuss the outcome of three alternative mechanisms of spool initiation: temporary pinning of the tip of the gliding
63 microtubule bundle at defective motors, simultaneous cross linking of multiple mic rotubules into a ring struct ure, and tip binding of motors. In a recent publication, Liu et al  proposed that spool formation results from bending of the tip of a microtubule or a microtubule bundle as a second microtubule wraps around it (due to the rotary motion of some microtubules)  We have reenacted this process with a variety of macroscopic tubula r structures and were unable to find any bending induced by the wrapping motion. Therefore we do not include this mechanism in our discussion. The newly developed understanding of the spool formation process can be utilized to optimize the active self asse mbly process, for example in order to favor the production of long bundles over spools. 4.3.2 Experimental Results After biotin/streptavidin covered microtubules adhere in random orientations to the kinesin coated surface, they are transported by the kines in motors at a velocity of 0.6 0.8 m s 1 and collide with and bind to each other forming extended linear bundles. ure 4 1 1 ). Several hundreds of these spools are imaged (in different fields of view and different experiments). While the imaging of the spools is typically conducted only after the formation of the spools has been completed in order to avoid interference of t he fluorescence excitation light with the assembly process, times lapse images are collected in some experiments.
64 From these time lapse data, it is apparent that the spool formation is often initiated by the pinning of the microtubule tip to the surface (F ig ures 4 11 B and 4 11 C ). These pinning events are well known  and typically attributed to a defective kinesin motor. In a previous study they occurred on average once every 750 m along the path of a microtubule  A second formation mechanism (Fig ure 4 11 D ) was observed only after being discovered i n simulations of the assembly process (simulations described below). Three (or more) microtubules approaching each other simultaneously can cross link in a triangular shape, which then relaxes into a circular shape over time. This process is favored in the first few seconds of the assembly process, when microtubule surface densities are highest. The length distribution of microtubules and microtubule bundles is determined by manually measuring the length of the structures observed in several experiments usi ng image analysis software. The size distribution of the spools which have formed after at least 20 min peaks at a circumference of about 10 m, and exponentially decays towards larger sizes (Fig ure 4 1 2B ). Spools with circumferences smaller than 6 m are not observed. The circumference is used as a metric here, since some of the larger spools are asymmetric, making the determination of a diameter more ambiguous. To obtain a sufficient number of spools for a statistical analysis, we have pooled the da ta from four separate experiments with streptavidin concentrations during assembly of 5 20 nM. The spool images of the four experiments were qualitatively the same. The appropriate choice for the microtubule persistence length is an interesting question, s ince the microtubule is a complex mechanical structure whose stiffness
65 depends on the experimental context [74,79 81] Pampaloni et al explained their observations of a length dependent persistence length by the limited longitudinal displacement between adjacent protofilaments enabled by bending lateral bonds  This limited displacement conveys high flexibility to a bending microtubule as long as the bending requires only small displacements between tubulins in adjacent protofilaments (applies to a large radius of curvature or a short microtubule). The micro tubule stiffens when bending increases as the lateral bonds cannot accommodate the displacement by bending alone. In our case, we believe that the large bending of the microtubule during spool formation results in a stiff response with a large persistence length of 5 mm  The persistence length of the microtubule also has to be distinguished from the smaller persistence length of the microtubule trajectory (0.1 mm) which results from small fluctuations of the microtubule tip only and is employed in the off lattice simulations described above. Under these assumptions for L p F, and w a minimal spool diameter of ~2 m (circumference of ~6 m) is obtained, in good agreement with previous observations [65,82] A novel simulation code  enables the simulation of hundreds of microtubules gliding on a surface and interacting with each other. In the simulation, microtubules represented as segment ed chains move on a triangular lattice with periodic boundary conditions and bind to each other at each collision (lattice size 2400 2400, segment size 80 nm, initially 1250 microtubules of 5 m length, results averaged over 200 runs with 3,000 seconds). Surprisingly, a random distribution of microtubules of average length and surface density similar to the microtubule population used in the experiment evolves into spools
66 with a size distribution which is very close to the experimentally observed size dis tribution (Fig ure 4 1 3 ), even if unlike in the experiments the microtubules are confined to movements on straight lines (unless they meet another microtubule and join it). However, spool formation is mostly complete within the first minute of simulated mic rotubule movement. Since in these simulations microtubules cannot turn by themselves, and since pinning events are not part of the simulation, the formation of spools in these simulations is entirely the result of simultaneous collisions between multiple m icrotubules. The obtained size distribution makes intuitive sense: it is unlikely that three or more microtubules meet at their tips and create a tiny spool; it is highly likely that they meet somewhere in their middles and create a spool with a circumfere nce of three times half the average microtubule length; it becomes increasingly unlikely that very long microtubules or microtubule bundles (which are rare) participate in the simultaneous collision to create large spools. 4.3.3 Discussion The experimenta lly determined distribution of spool circumferences (Fig ure 4 13) is in good agreement with previous observations [41,66,67,73,84,85] even though the m orphology of spools obtained by cross linking with streptavidin is somewhat different from the morphology of spools obtained by cross linking with streptavidin functionalized quantum dots [67,73] Specifically, the formation of broad spools with significantly different inner and outer diameters is less prominent here. Our off lattice simulations of the loop formation process (Fig ure 4 11 D Brownian bending) reveal that the selection bias towards smaller loops due to the limited length of the microtubule bundles is capable of producing surprisingly small spools
67 (circumferences of 100 20 0 m) compared to the average size of the loops in the trajectory. However, the experimentally observed spools are still smaller (circumference < 40 m), so that spool formation due to thermal fluctuations in the gliding direction cannot account for the ex perimental observations. Furthermore, the off lattice simulation is likely to overestimate the frequency of occurrence of small loops, because it does not consider the increasing microtubule stiffness due to prior bending suggested by Pampaloni et al  In contrast, the analytical modeling of a mechanism dependent on the mechanical work exerted by the motors, specifically the bending of the microtubule against the stationary t ip pinned by a defective motor (Fig ure 4 11 D pinning at defect), produces a picture in accord with the experimental observations. The absence of very small spools (<6 m circumference) is explained as a result of the inability of the motors holding the ben ding microtubule to provide sufficient attachment force. The declining frequency of larger spools is due to the declining prevalence of sufficiently long microtubule bundles. remain attached to the pinning motor during the entire spool formation process; an event which results in small spools. We have also observed release of microtubule tips from the pinning motor before a spool has fully formed, and only the subsequent collision of the bent front section with the center or tail section of the microtubule resulted in the formation of a larger spool. A more detailed modeling of the pinning process also has to account for the dynamic changes in the length distribution of microtubules and microtubule bundles within the first minute of the experiment. In
68 combination, t he increasing number of free parameters of such a model is not likely to lead to a proportional increase in our understanding. to spool formation (Fig ure 4 11 D simultaneo us sticking). Lattice simulations show that the cross linking of three or more filaments into closed, ring like structures which evolves into a circular spool can lead to a size distribution similar to the experimentally the third power of the microtubule surface density, which makes it most likely to occur before significant bundling of microtubules took place. Unfortunately, these first few seconds of the assembly pr ocess are difficult to image with our current experimental setup. The cross linking of two microtubules at their tips (Fig ure 4 11 D tip binding) can in principle give rise to a curved structure if the microtubule velocities are different, but spools of sma ll diameters are expected to result from this process. No distinct peak is observed in the experimental data, and the small probability of such tip binding events makes it unlikely that this is a major contribution to spool formation. frequent simultaneous collisions, favors the former, while low density, implying large distances between collision s and a high likelihood often countering a defective motor, favors the latter. However, both processes lead to similar spool size distributions, which make questions about the specific route less pressing.
69 The key insight is that spool formation is not act ivated by a Brownian ratchet type process, where rare and thermally activated bending events lead to spool formation. This result is important, because the fact that the formation of spools requires a significant amount of mechanical work to bend the micro tubules is in itself not proof that it could not be accomplished by thermally activated self assembly. The formation of biotin streptavidin bonds releases free energy, which could conceivably cause a microtubule or microtubule bundle to step by step (or bo nd by bond) wrap around itself once a spool has been initiated by thermal fluctuations. Spool formation by thermal fluctuations can indeed lead to spools much smaller than the persistence length of the microtubules, as our simulations have shown, but not s mall enough to explain the experimental observations. Instead, our investigation has determined that the primary Our second goal was to identify a mechanism to control the spool size. Manipulation of the system parameters, such as initial microtubule density and length as well as kinesin density, can be expected to modify the spool size distribution. However, since spool size and size distribution seems to be very similar indifferent ex periments under different conditions from different laboratories, the effects are likely to be small. This agrees with our analysis which shows that two very different formation mechanisms lead to very similar size distributions. The most productive approa ch to of defined sizes, for example in guiding channels 
70 The broad technology trends towards miniaturized device sand complex materials create the need for advances in assembly methods beyond therma lly activated self assembly (also known as chemistry when the components are molecules) and robotics. Active self assembly has the potential to make a contribution, and the kinesin/microtubule model system can be used to identify the principles underlying this approach. 4.4 Summary Biotinylated microtubules partially coated with streptavidin gliding on surface structures. We have presented a Cellular Automaton simulation techniqu e that models the kinetics of microtubule gliding and interactions. The analysis of the motion allows statistically meaningful data to be obtained which can be compared to experimental results, and aid in understanding the formation process of these nanost ructures. With this simulation tool, multiple scenarios of interest, including those without any current experimental realization, can be modeled quickly while varying the critical mechanistic parameters. Of course, the current simulations cannot capture a ll aspects of the experiments. The complex shape of spools is typically not observed experimentally, because spools tend to relax to a circular state to minimize their internal stresses. Experimental observations of a preference for a direction of rotation have also not be investigated, because the chirality of the microtubule lattice is not modeled here  However, a wealth of new ins ights into the assembly process has already been generated.
71 Figure 4 1 Image of a simulated nanospool with a circumference of 39.84 m, and composed of 449 microtubules all occupying the same lattice positions. Not all 449 microtubules and their posi tions are reflected in the simulation in order that attention can be focused on the formation of the nanostructures. Regions of overlap from the microtubules are expressed above in the multitude of colors within the nanospool.
72 Figure 4 2 Initial, intermediate and final state of a representative simulation: Initially 200 microtubule filaments are randomly arranged on an 80 m x 80 m hexagonal lattice. At time step 350 (35 seconds) some nanowire and nanospool structures have formed. At time step 50,000 (5,000 seconds) the simulation is complete and only nanospool structures remain
73 Figure 4 3 Formation of a simulated nanospool. At time step 10 (1 second), A nanowire is approaching intersection with itself at a body subunit. At time ste p 12 (1.2 seconds) intersection takes place. Upon intersection, the head subunit corresponds to the point of intersection where the microtubule merges to begin to form a nanospool. There is an overlap in lattice positions.
74 Figure 4 4 Confirmation of the turn probability values calculated from microtubule movement within the simulation with the initial condition parameters set by the users. Within this example, the microtubule is set to take on 60 o turns 4.5% of the on the hexagonal lat tice. (All microtubules are restricted in moving in the backwards direction on itself. ) The sum of the total movements ( forward direction 60 o turns, and 120 o turns ) should equal the product of the number of microtubules and time steps. The percentage of t he number of sixty and one hundred twenty degree turns should equal the initial turn probability percentage dictated by the user
75 Figure 4 5 Analysis of the output from 200 simulation runs on a grid of size 200 m x 200 m with 1250 microtubule chain s, a uniform 5 m microtubule length, and a turn probability of 0.045%. The number of nanospools formed per run fluctuates and can be fitted with a Poisson distribution. The distribution of spool circumferences peaks at 12 m. A second peak near a circumfe rence of 200 m is due to the periodic boundary conditions and stems from spools spanning the whole simulation area. The radial distribution function describes the distribution of spool locations on the surface. Radial distribution functions are calculated for each of the 200 runs.
76 Figure 4 6 The effect on the simulation output as the MT density is varied. The simulation output includes: (a) number of spools generated; (b) spool circumference; (c) radial distribution function is analyzed. The simulation was performed on a 200 m x 200 m surface for 10,000 seconds.
77 Figure 4 7. The effect on the simulation output as the MT length is varied. The simulation output includes: (a) number of spools generated; (b) spool circumference; (c) radia l distribution function is analyzed. The simulation was performed on a 200 m x 200 m surface for 10,000 seconds.
78 Figure 4 8. The effect on the simulation output as the MT turn probability is varied. The simulation output includes: (a) number of spool s generated; (b) spool circumference; (c) radial distribution function is analyzed. The simulation was performed on a 200 m x 200 m surface for 10,000 seconds.
79 Figure 4 9 Analysis of the number of nanospools generated versus time from one simulation trial on surface with 200 m x 200 m dimensions, 1250 microtubule chains, a uniform 5 m microtubule length, a nd a turn probability of 0.045% : (a) the number of spools generated in the first 10,000 seconds ; (b) the number of spools ge nerated within the first 200 seconds. The number of spools formed was fit with a logistic fit.
80 Figure 4 1 0. Description of previous studies describing interactions proposed to initiate nanospool structures ( a ) Microtubules are polymerized from tu bulin dimers functionalized with biotin linkers. Partial coverage of these linkers with tetravalent streptavidin enables cross linking of microtubules. ( b ) Kinesin using ATP as a source of chemical energy. ( c ) Collisions between sticky microtubules lead to the formation of elongated bundles and finally spools. ( d ) Initiation of spool formation from microtubule bundles can potentially result from thermally activated fluctuations i n the direction of the microtubule movement (Brownian bending), pinning at defective kinesin motors, simultaneous aggregation into a ring like structure, or tip binding of microtubules moving at different velocities. Reproduced from Ref erence 
81 Fig ure 4 1 1. Experimental images of nanospool formations ( a ) Aggregation of microtubules leads to the formation of spools of different sizes. ( b and c ) Pinning of the leading tip of a microtubule bundle initiates spool formation. ( d ) Simultaneous collisions of three (or more) microtubule bundles form triangular structures relaxing into circular spools. Reproduced from Ref erence 
82 Fig ure 4 1 2. Experimental size distribution of spool circumference was measured for 607 spools. Reproduced from Ref e rence 
83 Fig ure 4 1 3. Comparison amongst experimental and simulation nanospool circumference length distributions Initial ( a ) and final (b) snapshots of a run of the lattice simulation of microtubule s pools formation. (c) Comparison between spool circumference distributions of spools formed in the lattice simulation (red hashed bars) and experimentally observed distribution of spool circumferences (grey bars, normalized). Reproduced from Ref erence 
84 CHA PTER 5 POLYTHIOPHENE THIN F ILM GROWTH VIA SURFA CE MODIFICATIONS 5 1 Thin Films Collision induced reactions are commonly used to treat material surfaces. They involve an incident particle and a target, where the incident particle is composed of functional materials of interest and the target is thiophene. The particles are typically a ccelerated at high velocities with incident energies that range from fractions of an eV to 10 MeV [87 89] Particles accelerated at hyper thermal energies typically come from sources such as plasma  laser beams  or ion beams  A variety of interactions are possible, including the chemical modification of the surface and the deposition of thin films [93 96] The properties of the thin films are engineered by controlling how the parti cles approach and interact with the substrate. 5 2 Growth of Organic Films Organic thin films are used in a range of applications such as light emitting diodes, field effect transistors, photovoltaics, sensor films, recording materials and rechargeable batteries  The films are commonly grown from the vapor phase by thermal and physical deposition. Thermal deposition methods include: plasma, laser beam, or ion beam methods [71,97] Physical deposition involves solution based methods such as spin coating, casting, and printing [98 100] Direct thermal deposition uses low molecular weight oligomers or polymers which degrade and produc e gaseous species These species deposit on the substrate surface yielding the thin film. In other physical so lution based methods, thin films are generated through the self assembly of molecules [101,102]
85 Both direct thermal deposition and solution based physical deposition methods grow thin films but each method has its own set of limitatio ns. For example, small molecules from the degradation of polymers and oligomers condense from the gaseous state into the chambers of the thermal deposition apparatus due to their low molecular weight. In addition, the thickness of thin films is difficult t o manage in solution based physical deposition methods due to a lack in control of ordering of the molecules. The thickness of thin films usually ranges from a monolayer at a few nanometers; however, the thin film thickness produced from solution based met hods surpasses a few monolayers. The maximum thickness of a thin film is several micrometers  5 3 Organic Material Surface Modification As deposited material comes into contact with the target surface, interactions take place that chemically modify the surface. Modifications are primarily made using plasma and ion beam deposition [103 105 ] Plasma, an ionized gas with a minimum of one dissociated ion from an atom, molecul e, or cluster, is generated from the collisions of accelerated particles [94,106] Accelerated particles react with the ta rg et al lowing components to deposit and modify the target. Plasma treatment offers advantages to thin films where modifications can be confined to the depth of several hundred angstroms in the surface layer without modifying the bulk properties of the poly mer. Additional advantages include the fact that alterations made in the surface are fairly uniform and changes can be made to the surfaces of all polymers, regardless of their structures and chemical reactivity  However, there are disadvantages to plasma deposition where specie identification becomes a challenge due to the complex environment in which numerous association and dissociation reactions take place. Additional disadvantages include increased cost of operation due to the need to operate
86 under vacuum conditions and the lack of uniform settings that could be used for treatment in all ma terial systems  Ion beam deposition performs the same function as the plasma in the deposition of thin films onto the target but a different mechanism is used. High energy ions are created under vacuum by an ion source and accelerated by a voltage to form a beam  The ions collide with the atoms of the polymeri c material producing vacancies in the outer layer of the polymer Consequently, the ions become implanted in the vacancy locations. Consequently, an outer layer forms yielding a thin film  T he energy of the ion beam can be increased past the vaporization point, to permit the removal of the material to later redeposit and form a thin film The microstructure of the deposited film can have novel properties such as reduced grain size, the presence of metastable phase s and improved mec hanical properties  Ions also transfer energy to a material as interactions between the ions and materials occur. Ion beams at low energies do not penetrate deeply into polymer surfaces, whereas at high energies the ions can penetrate up to several micrometers into a treated polymer  For example, an ion beam with energies in the range of a few hundred kilovolt range penetrate a treated material to the depth of 1 10 m  However, there are additional factors that affect the penetration depth of ions in addition to their energies. T hese include the ion mass and the nature of the polymeric material (molecular weight, crystallinity, and cross linking). In general, ions with the highest energy and lightest mass possess the highest penetration depth  Ion beams off er the advantage of low processing temperatures, prevention of heat damage, and yielding localized polymer surface modifications. There is additional
87 control in the ion specie composition of the resulting beam, where desired ions used are isolated and sort ed by mass selection for surface modification and deposition treatments in a material  On the other hand, ion beams are very energetic and can modify the polymer surface properties, leading to excessive surface modifications, unnecessary chemical reactions, and possible surface damage  5 4 Surface Polymerization by Ion A ssisted Deposition The surface polymerization by ion assisted deposition (SPIAD) method is an alternative experimental technique used to combat the limitati ons from thermal and physical deposition in the deposition of thin films. Most impor tantly, SPIAD can manipulate thin film s grown to yield desired properties. Hyper thermal polyatomic ions, such as C 4 H 4 S + C 3 H 3 + CHS + C 2 H 2 S + C 3 H 5 + and neutrals are simultaneously deposited in vacuum to form the thin films [13,112] SPIAD offers flexibility in the use polyatomic ions because they can be mass selected or non mass selected. Mass selected ions are used in studies of physical forces, wh ere ions are controlled. Non mass selected ions are used in prototype manufacturing processes. SPIAD offers an advantage in producing a multitude of thin films that can be engineered by varying features such as the ion energy, ion structure, ion/neutral ra tio, neutral structure, and substrate temperature  Properties of the thin film properties such as the morphology, electronic structure, film thickness, and properties of the target material can be optimized [14,15,112] SPIAD allows fine tuning of the thin film opti cal band gaps, alterations in the optoelectronics of the thin film through chemical changes, and morphological changes [13,15] This dissertation will examine collision induced reactions in the gas phase to focus on the surface polymerization of polythiophene thin films by ion assisted deposition.
88 Polythiophene, a conductive polymer, is attracting much interest due to its ability t o remain stable at high temperatures and because its molecular structure can be readily transfigured to achieve desired electronic properties  Chemical reactions of thiophene with organic molecules are of interest to chemically modify thermally deposited coatings or thin films of conductive polymers. However thermally deposi ted polymeric thin films degrade in response to mechanical deformation and extreme thermal fluctuations [8 10] Chemical modification is one way to stabilize these films without substantially altering their electronic and molecular structure. Modifications have been ma de to thiophene with SPIAD by Hanley and co workers [11 15] In this approach, thiophene ion deposition is combined with simultaneous dosing of terthiophene vapor to achieve the app ropriate level of modification while maintaining performance. In particular, beams of mass selected thiophene ions at 60 nA are produced by electron impact, which causes chemical reactions to take place. This method produces polythiophene films with desire d optoelectronic properties that largely  The SPIAD approach provides control over polymerization and thermal film growth through sputtering, bond breakage, energe tic mixing, and other processes [15,114] The initial impact of SPIAD allows for the polyatomic ions to interact only with the top few nanometers of the surface. As a result, there is tremendous control of the nanoscale morphology [114,115] Importantly, SPIAD results in selective polymerization to stabilize the thin film, where the polythiophene films produced display fluorescence in the UV/vis region, an important property of conducting polymers 
89 Ab initio methods will be used to investig ate reactive systems involving thiophene and small hydrocarbon radicals. The results will provide insights into the chemical processes that occur in the chemical modification processes described above. In particular, a series of calculations within the GAU SSIAN03 program  are used to calculate the thermochemi stry of reactions from the enthalpy of formation, and to identify the corresponding transition states. In addition, we consider the influence of the level of theory used in the results by evaluating the performance of various hybrid functional methods.
90 CHAPTER 6 HYBRID FUNCTIONAL ME THODS 6.1 Basic Quantum Mechanics to the Schrdinger Equation 6.1.1 Time In dependent Schrdinger E quation The electronic structure of a material is determined by calculating the energy associated with the collection of atoms that make up the material, their electrons, and its structure. Solving the time independent Schrdinger equation (Equation 6 1) allows the lowest energy state, the ground energy, to be determined, where H is the Hamiltonian operator, is the electr onic wave function of the spatial coordinates of all electrons, an d E is the ground state energy. (6 1) The wave function describes the properties of many electrons; the wave function of a single atom cannot be identified without considering the interactions of all the electrons associated with the atom with each other. The standard interpretation is that the product of (where the superscript denotes the complex conjugate) represents the spatial distribution of the electron probabilit y density 6.1.2 Born O ppenheimer A pproximation The Born Oppenheimer Approximation  can be applied to the time independent Schrdinger equation in order to solve for the wave function. The mass of a nucleon (proton or neutron) is ~1800 times that of an electron. The approximation allows for the separation of the electronic and nuclear co mponents. The positions of the nucleus are assumed fixed  and their interactions and kinetic energy is not calcul ated. Instead, the motion of the electrons are accounted for  and the ground state energy of the electrons are expressed as a function of the atomic positions 
91 In the time independent Schrdinger equation, the Hamiltonian is now represented as three components: the kinetic energy of each electron, the interaction between each electron and the collection of atomic nuclei, and the interaction energy among the electrons respectively, as given in Equation 6 2  The mass of an electron is represented as me; i and j refer to the electrons, A represents the nuclei, Z represents the charge of the nucleus is the radial position of an electron, and N is the total number of electrons. = ( 6 2) The potential energy surface diagram can be extracted as the nuclei positions are varied in small steps and the Schrdinger equation is solved for the ground state energy. Without the Born Oppenheimer ap proximation ther e would be a lack in the concept s of the potential energy surface, equilibrium, and transition state geometries  6.2 Approaches to Approximate Solutions in the Schrdinger Equation 6.2.1 Hartree Product Single electron approximations are used in the Hartree product to calculate the wave function for in an N electron system. The total wave function is represented as eigen function s of single electron spin orbitals in Equation 6 3 which includes the spatial coordinates of electron positions, x 1 where i and j denote the specific electron and its spin state for all N electrons (6 3) The Hartree product is not a good approximation of the total wave function because it is not antisymmetric as required for a Fermion 
92 6.2.2 Hartree Fock (HF) Approximation The Hartree Fock HF, approximation incorporates antisymmetry among interchange of electrons that the Hartree Product lacks. In the Hartree Fock method the wave function is approximated as single electron functions where the exact solution of the wave function is determined from the Schrdinger equation  More specifically, the complete wave function is approximated as a single Slater determinant. A Slater determinant is used in Equation 6 4  to obtain a better approximation for the wa velength. (6 4) The determinant of a matrix is solved for the single electron wave functions as shown in Equation 6 4 for two electrons. In the equation there is a normalization factor of and a built in electron exchange implicitly where a sign change occurs each time two electrons are exchanged  All of the single particle wave functions in the N electron system are combined according to the antisymmetry rule, which is the Pauli Exclusion Principle  In HF, the total electron wave function is a Slater determinant of the spin orbitals from the N single electron equations. The spin orbitals are approximated in Equation 6 5 using expansion coefficients, j,i and a finite set of functions, i (x), also known as a basis set. (6 5) An infinite basis set would converge to the exact wave function but at an infinite computational expense. Therefore, in practice, a finite basis is used. A larger basis approximates the true wave function better but it is computationally costly; a smaller
93 basis set is more computationally efficient but worse at approximatin g the true wave function. The Hamiltonian in the HF Equation (6 6) is composed of three terms: term one represents the Coulombic interactions amongst nuclei, where A and B denote the nucl ei (this term does not depend on the electrons); term two is the sum of the kinetic energy of the electron and the potential energy of the Coulomb attraction between the nuclei and electron n This term is defined as ; ter m three is the Coulombic repulsion interactions amongst e lectron pairs. (6 6) Figure 6 1 highlights the iterative procedure used to solve for the wave function in the Hartree Fock Approximation, while Equation (6 7) combines all of the components within the Hartree Fock Approximation to solve for the ground state energy in the Schrdinger e quation, where is the potential energy that accounts for repulsion between nuclei (Equation 6 8) is the kinetic energy for non interacting electrons (Equation 6 9), is the potential energy of the Coulombic attraction between electrons and nuclei (Equation 6 9) and is the potential energy of the Coulo mbic interactions of electron electron interactio ns (Equation 6 10) (6 7) ( 6 8 ) = = ( 6 9 ) 1 2 3 (6 10) ( 6 9 ) ( 6 8)
94 ( 6 10 ) 6.2.3 Density Functional Theorem (DFT) Methods: Hohenberg Kohn and Kohn Sham Density Functional Theory, DFT, is an approach used to understand the ground state properties of matter at the level of individual atoms  Three spatial variables for each electron total to 3N variables to be taken into account for a system  The summation of these interactions causes the calculation to be a many bodied problem. DFT is centered on its use of functionals, a function of a function, of the electron density. The reduction in computation time that this provides offers advantages in electronic calculations; however, the system size of a material is limited to a several hundred atoms. DFT is able to probe a large a variety of details, including structures, vibrational frequencies, atomization energies, ionization energies, electric and magnetic properties, reaction paths, and enthalpies of formation. 184.108.40.206 Hohenberg Kohn theorems Pierre Hohenberg and Walter Kohn provided a mechanism to solve th e Schrdinger e quation for the electron density instead of the wave function. There are as many electron wave functions as there are atoms; however, there is only one atomic density; both are functions of position This is done to minimize the number of variables to solve for in the Schrdinger e quation In solving for the wave function there are 3N variables, and the electron density accounts for 3 spatial coordinates. As a result of building the approach around t he electron density, the ground state energy can be calculated, thereby providing an approximate solution to the Schrdinger equation. The ground state energy is represented as a functional of the electron density, where E Ground State = E [n(r)]. The ground state energy is determined from the electron density due to a
95 1:1 mapping between the ground state energy and the electron density  The electron density that minimizes the energy of the overall functional defined in Equations 6 11 and 6 12 is the true electron density that solves for the full solution of the Schrdinger e quation  (6 11 ) (6 12 ) The ground state energy is defined to have known and exchange correlation components in Equation 6 8. The known energy component in Equation 6 9 accounts for the kinetic energy of electrons, Coulombic interactions amongst the electrons and atomic nuclei, Coulombic interactions between electron pairs, and Coulombic in teractions amongst nuclei pairs. The exchange correlation functional, in E quation 6 8 is not uniquely defined but it is recognized to include all quantum mechanical effects that are not included in the known energy component  6 .2. 3 .2 Kohn Sham e quation The Kohn Sham equation ( 6 1 3 ) provides a mechanism to find an approximation to the correct electron density with the purpose of minimizing the total energy of the functional The single electron wave functions depend only on three spatial variables, (r) The components of the Hamiltonian in the Kohn Sham equation respectively include : the kinetic energy of electrons (Equation 6 1 3 ) ; a C oulombic interaction amongst an electron and a collection of atomic nuclei (Equation 6 1 4 ); the Hartree potential (Equation 6 1 5 ) represents a Coulomb repulsion between an electron and the total electron density by all electrons in the system. Part of includes a self interaction contribution because the electron contributes to the total electron
96 density; the exchange correlation term that represents interactions amongst the electrons (Equation 6 1 6 ) The exchange correlation term is the functional derivative of the exchange correlation energy. It is used as a correction for self interactions amongst el ectrons, and provides correlation and exchange to single electron equations. An iterative process outlined in Figure 6 2, involving the electron density, the Hartree potential, and the single electron wave functions, is used to solve the Kohn Sham equation s  The ground state is determined by m inimizing the energy of the energy functional in Equation 6 11 to find a self consistent solution to single electron equations  (r) = (6 1 3 ) V(r) = (6 1 4 ) where n(r) = (6 1 5 ) (6 1 6 ) 6.2.4 Relative Advantages and Disadvantages of HF vs DFT methods Hartree Fock methods offer the advantage over DFT in providing the exact solution of the Schrdinger equation to a simplified problem. However, a major weakness t o the Hartree Fock Approximation is that it neglects electron correlation, which is physically incorrect and can lead to large deviations from experimental results  Post Hartree Fock methods have been developed to include electron correlation to the multi electron wave function. For example, M ller Plesset perturbation theory treats correlation as a perturbation 
97 An uncertainty is always present in DFT calculations since exact solutions are not known  DFT is based on the electron density, which is a property that has a clear physical realization, as compared to the electron wave functions whose physical interpretation is less obvious. The wave function in the Hartree Fock method becomes more complicated mathematically as the number of electrons increases. In DFT, the density depends only on three variables that are based on the spatial x y z coo rdinates of the individual electron scaling N 3 instead of N 4 from the wave function, where N represents the total number of electrons in a system  The fewer variables allows DFT to be computationally efficient. There are a number of different DFT methods with different forms for the function al. Some functionals work well for some systems; while others work best for other systems. For example, the Becke 3 Parameter (Exchange) Lee, Yang and Parr correlation (B3LYP) functional performs in an appropriate manner for transition metal applications but not for organic compounds [122 125] DFT also limits the number of atoms to examine in a system to a few hundred due to computation expense. DFT also generates inaccurate results when there are weak wan der Waals attractions between atoms and molecules 6. 3 Hierarchy of Approximations The exchange correlation functional must be defined in the Kohn Sham equations in order to solve for the electron density in the Schrdinger equation; however, defining this term is difficult. Different approximations of the e xchange correlation functional are used in the Kohn Sham equations. More specifically, there is a hierarchy of functionals which include more detailed information. The hierarchy from the least physical information added toward the method that can comes clo sest to an exact solution to the Schrdinger equation includes Linear Density Approximation (LDA), Generalized
98 Gradient Approximation (GGA), Meta Generalized Gradient Approximation (Meta GGA), and Hyper Generalized Gradient Approximation (Hyper GGA). Each approximation in the hierarchy has several functionals in each stage  LDA approximates the exchange correlation from the local electron density. The electron density is assumed to be constant and the exchange correlation potential is defined in Equation 6 1 7 where the exchange correlation potential for a spatially uniform electron gas is the same density as the local electron density. The GGA functional incorporates more information than the LDA functional ; the electron density is slowly varied and the correlation function accounts for the local electron density and the gradient in the electron density in Equation 6 1 8 The meta GGA functional differ s from the GGA functional with the addition the Laplacian of the electron density in Equation 6 1 9 The kinetic energy of the Kohn Sham orbitals as shown in Equation 6 20 contains the same information as the Laplacian of the electron density, and can be used as a substation in some meta GGA functions. Hyper GGA functionals are composed of a mixture of the exact exchange derived from the exact exch ange energy density in Equation 6 21 and a GGA exchange functional. The exchange functional is nonlocal and evaluated over all spatial locations. [n(r)] (6 1 7 ) [n(r), n(r)] ( 6 1 8 ) [n(r), n(r), 2 n(r),] (6 1 9 ) (6 20 ) = (6 21 )
99 6 .3. 1 Hybrid Functionals F unctional s that incorporate DFT and HF are r eferred to as Hybrid functional s There are two categories for Hybrid functionals that will be the main focus of discussion in this chapter: Hybrid GGA and Hybrid meta GGA functionals. Hybrid GGA functionals comb ine an exact change with a GGA functional, while the Hybrid meta GGA functional incorporates exact change with a meta GGA functional. The Hybrid GGA functionals of interest will be the Becke 3 parameter Lee Yang Parr (B3LYP) Hybrid functional and the Becke 's 1998 revisions to B97 (B98) functional, and the meta GGA functional is the dependent hybrid (BMK) functional. 6.3.2 Becke 3 parameter Lee Yang Parr (B3LYP) H ybrid Functional The Becke 3 parameter Lee Yang Parr (B3LYP) Hybrid functi onal is the most common ly used hybrid functional due to its flexibility for use in several environments. In the B3LYP functional, the exchange correlation energy is composed of the four terms in Equation 6 22  The first term is the local spin density (LSDA) factor  of the exchange energy density, which is the Hartree Fock exchange energy of an inhomogeneous many electron system of a uniform electron gas (Equation 6 2 3 ) The signifies the electron spin state  The LSDA factor usually underestimates the exchange energies for atomic and molecular systems by approximately 10%; The second term replaces some o f the electron gas exchange with the exact exchange, where the represents the exact exchange energy ( Equation 6 24) 1988 gradient,  accounts for correction to the LSDA factor in the third term (Equation 6 25) ; The fourth term, is the 1991 Perdew and Wang gradient correction for correlation, where = 0.004235, = is the local polarization of
100 the electrons, and n is the density of electrons (Equation 6 2 6 )  Physical information is embedded into the B3LYP functional with semi empirical coefficients, which are fit to experimental data. The values of the coefficients are respectively listed as: 0.2, 0.72, and 0.81  = (6 22 ) (6 2 3 ) (6 2 4 ) = (6 2 5 ) = (6 2 6 ) dependent (BMK) Hybrid Functional The Boese dependent hybrid functional (BMK) was developed to examine reaction mechanisms  The BMK exchange correlation energy has two exchange terms, (Equation 6 28) and (Equation 6 29 ), a correlation term, (Equation 6 3 0 ) and a term that includes the product of the Hartree Fock exchange energy  and a mixing coefficient, (Equation 6 31 ). In the BMK functional, the zeroth order exchange coefficient (Equation 6 32) varies from 1.1 ( at 0% HF exchange ) to about 0.4 ( at 50% HF exchange ) and the coefficient (Equation 6 33) changes from near zero (0.001 ) at 0% HF exchange to 0.38 at 50% HF exchange. Both coefficients contribute most to the overall exchange correlation energies. Details regarding the kinetics of a reaction are extracted through the kinetic en ergy density, ( 6 23 ) ( 6 24 ) ( 6 25 ) ( 6 26 )
101 which is defined Equation 6 3 4 is used and contributes to the exchange correlation energy. Equations 6 35 through 6 40 define variables that are components of E quations 6 28 through 6 34. (6 2 7 ) (6 28) (6 29) (6 30) = (6 31) (6 32 ) = (6 33) = (6 34) (6 3 5 ) = (6 3 6 ) = (6 3 7 ) (6 3 8 ) = (6 3 9 ) (6 40 ) ( 6 28 ) ( 6 29 ) ( 6 30 ) ( 6 31 )
102 onal (B98) is composed of exchange and correlation terms (Equation 6 41 ). The exchange energy, is approximated from the Perdew [ 131] and Becke  functionals in Equation 6 42 The local spin density (LSDA) factor (Equat ion 6 23 )  spin exchange energy density of a uniform electron gas with a spin density equal to The gradient corrects the LSDA factor through the dimensionle ss density gradient, The correlation term in B98 is composed of two terms listed in Equation 6 43 The first term, (Equation 6 44 ) represents antiparallel spins, and the second term, (Equati on 6 45 ) represents parallel spins from uniform electron gas in the correlation energy. Both terms are calculated from the Perdew and Wang parameterization  in equations 6 46 and 6 47 Two gradient dimensions, and are dependent on the dimensionless variables and (Equation 6 48) (6 41 ) ( ) r (6 42 ) (6 43 ) = (6 44 ) = (6 45 ) = (6 46) ( 6 42 ) ( 6 43 ) ( 6 44 ) ( 6 45 )
103 = (6 47 ) (6 48) 6.4 Calculation Details All of the calculations discussed in this chapter were performed with the Gaussian 03 software package, a modeling program that examine s the electronic structure of materials  There are user defined parameters, such as the functional and basis set, withi n Gaussian to aid in examining systems. The basis set parameter is critical because it approximates the wave function. Expansion of a basis set increases accuracy in calculations (due to a more accurate wave function approximation) and boosts computational expense. Additional time is needed for the basis set to converge to a solution. The basis set of 6 311+G (3df, 2pd), a relatively large basis set, is selected for use in the enthalpy of formation and transition state calculations. The 6 311+G specifies th e use of the 6 311G basis set for first row atoms, the McLean Chandler basis set [133,134] for second row atoms, and the supplementation of diffuse functions  The second portion (3df, 2pd) represents 3 sets of d functions and one set of f functions on heavy atoms, and supplemented by 2 sets of p functions on hydrogen atoms. The enthalpies of formation are calculated using a range of hybrid functional theories, including B3LYP  ; BMK  ; and B98  In addition, the G2  and G3  theories are used to calculate enthalpies of formations at 0 and 298K in E quations 6 49 and 6 50 M represents the molecule of interest, x is the number of constituent atoms, is the atomization energy (Equation 6 51) is the total electronic energy of the molecule, and is the zero point vibration energy of the
104 molecule. The and values are calculated in frequency calcul ations in Gaussian. The values are extracted from the thermochemistry. (6 49 ) (6 50) (6 51) The G2 method provides a n adoption to a higher level of correction with a better fit to experimental atomization energies. The G3 method improves accuracy further through a sequence of single point energy calculations using different basis sets, a new formulation of the higher level correction, spin orbit correction for atoms, and correction for core correlation Lastly, the Complete Basis Set Method, CBS QB3  is also used All theories, with the exce ption of CBS QB3, use the 6 311 + G ( 3df, 2p) basis set The Gaussian, G2 and G3, and Complete Basis Set, CBS QB3, met hods are known to produce highly accurate values for thermochemical parameters [138,139,141 143] Therefore, the G3 (for smaller systems) and CBS QB3 (for larger systems) metho ds are used as standards for systems when experimental data is not available. T ransition states are identified for each reaction The internal reaction coordinate (IRC) calculation is performed to verify the existence of the transition state by examining t he reaction path that connects the desired reactants to the expected products The reactivity of the hydrogen atoms at the ( the hydrogen bonded to the carbon atom closest to the sulfur atom in the aromatic ring) and (the hydrogen bonded to the second carbon atom from the sulfur atom in the aromatic ring) sites of the thiophene
105 molecule are considered. The B3LYP functional is selected for use in all the transition state calculations. The basis set of 6 311+G(3df,2pd), a relatively large basis set, is us ed for all enthalpy of formation and transition state calculations (with the exception of the CBS QB3 enthalpy of formation calculation s). Sections 6 .4.1 6 .4.2 and 6.4.3 briefly discuss the G 2, G3, and CBS QB3 methods. 6 4.1 Gaussian 2 ( G2 ) Method The G2 method is composed of seven calculations to obtain several pieces of information on the electronic structure of a material  : (1) a Mller Plesset ener gy correction truncated at the second order (MP2), optimization calculation to obtain the molecular geometry; (2) a quadratic configuration interaction calculation (QCISD(T)) with single, double, and triple excitations contributions to solve the nonrelavistic Schrdinger equation within the Born Oppenheimer approximation This cal culation is the highest level of theory, where interactions of different electronic configurations are examined requiring long computational time (and expense); (3) a Mller Plesset energy correction, truncated at the fourth order (MP4) calculation to asse ss the effect of polarization functions; (4) a MP4 calculation to evaluate the effect of diffuse functions; (5) a MP2, optimization with a large 6 311+G(3df,2p) basis set; (6) a HF geometry optimization calculation; (7) a frequency calculation to obtain de tails on the thermochemistry of a material by obtaining the zero point vibrational energy (ZPVE) The ZPVE is scaled by 0.8929 An empirical correction is added to account for factors not considered. This is called the higher level correction (HLC) and is given by 0.00481 x (number of valence electrons 0.00019 x (number of unpaired valence electrons).
106 6 4.2 Gaussian 3 ( G3 ) Method The G3 method is very similar to G2 but additional details are added. The 6 311G basis set in calculations 3 5 is replaced by the smaller 6 31G basis. The final MP2 calculation in calculation 5 use s a larger basis set, which is known as G3large All electrons are correlat ed unlike the valence electrons which are correlated in G2 As a result, core correlation contributions are ad ded to the final energy. The HLC has a same form as in G2 but different empirical parameters are used in the calculation  6.4.3 Complete Basis Set ( CBS QB3 ) Meth od The Complete Basis Set (CBS QB3) method is similar to the G2 and G3 methods; however the CBS QB3 method contain s a MP2 extrapolation in one step.
107 Figure 6 1. A flow chart of the iterative process to find a solution in the Hartree Fock Approximation The solution to the HF Approximation involves solving for the orbitals and ultimately the wave function in the Schrdinger equation. A n initial guess is made on the orbit al s by specifying an expansion coefficient; the electron density is estimated and then the singe electron equations are solved from the electron density to solve the spin orbitals If the spin orbitals found are consistent with the initial guess orbitals t hen these are the solutions If not, then a new estimate is made to the electron density and the process starts again until a solution is found A basis set is used to expand the orbital. Using a large basis set increases accuracy in the calculation and th e amount of effort to find a solution.
108 Figure 6 2. A f low chart of the iterative process to find a solution to the Kohn Sham e quations A solution to the Kohn Sham equation involves solving for the ground state electron density and ultimately the total ground state energy in the Schrdinger equation. A n initial guess is made for the electron density The trial electron density is then use d in the Kohn Sham equations to solve for the single particle wave functions, (r). The electron density is then calculated from the Kohn Sham single particle wave functions using the equation = 2* If the electron density from the Kohn Sham single particle wave functions is the same as the trial electron density used, then the electron density is the true ground state electron density that can be used to calculate the total ground state energy in the Schrdinger equation. I f the electron densities are different, then the trial electron density is updated and the process begins again by solving the Kohn Sham equations with the updated electron density.
109 CHAPTER 7 SURFACE MODIFICATION REACTION RESULTS 7.1 Enthalpies of Formation The enthalpies of formation for gas phase compounds are used to examine the structures in reactive systems. In particular, interactions involving thiophene with organic molecules are observed to determine the amount of energy loss and absorbed. T he enthalpies of formation are calculated from the thermochemistry of different levels of ab initio theories in order to determine which method most closely agree with the experimental values. The Gaussian, G2 and G3, and Complete Basis Set, CBS QB3, metho ds are known to produce highly accurate values for thermochemical parameters [138,139,141 143] Therefore, the G3 (for s maller systems) and CBS QB3 (for larger systems) methods are used as standards for systems where experimental data is not available. The calculated enthalpies for the molecules in reactions 1 4 (shown in Figure 7 1) are expressed in Table 7 1. For the mole cules considered here, the BMK and B98 approaches differ from the B3LYP functional by similar amounts. Of the eight systems tested, the B98 and BMK functionals both outperform the B3LYP functional in three cases each. The calculated values using B3LYP, BMK and B98 are in agreement with published values [139,144,145] The B98 f unctional performs best for large molecule systems, including those involving the thiophene radical and dithiophene molecules. In particular, the error ranges from 0.12 1.83% with the G3 method as the standard. In contrast, the BMK functional yields opti mal results predominantly for small molecule systems, including those involving CH2, CH3, and thiophene. In these cases, the error ranges from 0.36 2.29% in comparison with experimental data. The B3LYP functional
110 excels for the C 2 H and C 2 H 2 molecules, wi th errors ranging from 2.29% to 3.42%. Within the limited number of molecules considered here, this comparison makes clear that the enthalpies of formations calculated with BMK and B98 functionals yield reasonably accurate results. 7.2 Reaction Pathways The predicted reactive pathways for reactions 1 and 2 are shown in Figures 7 2 and 7 3, respectively. A C 2 H and CH 2 molecule, respectively, is added to thiophene resulting in the extraction of a hydrogen atom from the thiophene. In reaction 1, the C H bond from the thiophene molecule elongates and fully detaches. In reaction 2, the C H bond is elongated and eventually disconnects. The transferring hydrogen atom and the CH 2 molecule are nearly coplanar. In reaction 3, two thiophene radicals interact in a bar rierless reaction to polymerize and produce dithiophene, C 8 H 6 S 2 A thiophene molecule reacts with a thiophene radical in reaction 4 yielding a barrierless reaction and produces C 8 H 7 S 2 In comparing the reaction paths, it is conclude d that an increase in th e reactivity of the reactant molecules, such as thiophene + C 2 H, yields a high er activation energy. The site of the thiophene molecule is predicted to be the preferred location for the reactions due to the lower energy from the barrier heights for the tr ansition states and products. The hydrogen atom at the site is less reactive than the hydrogen at the site for two reasons. First, the electronegativity from the sulfur atom produces a slightly positive charge on the carbon atom. This site prefers no t to break the carbon hydrogen bond to abstract the hydrogen atom. However, if this does occur additional energy will be needed. Second, electron repulsion will occur at the site amongst the
111 products formed between the thiophene radical and the slightly negative charged sulfur. The site provides more distance between electron rich species. The IRC method  is used to verify and connect the reaction paths with the reactants and product. In reaction 1, the potential energy surface indicates an attraction between the C 2 H and thiophene reactants, which leads to the formation of a fully coordinated intermediate reactant complex (Figure 7 2) Therefore, there is a decrease in energy to a minimum. The reactant complex overcomes energy barriers of 57 kcal/mol (2.5 eV) for the site reaction and 45 kcal/mol (1.9 eV) for the site reaction. The imaginary frequencies indicating the transition states were 477.6 cm 1 ( site) and 198.9 cm 1 ( site). Reaction 1 yields isolated products of C 2 H 2 and a thiophene radical. Herbst et al found a large energy barrier in a similar reaction involving CH 2 and C 2 H 2  An intermediate complex forms among the reactants and a 38.8 kcal/mol energy barrier has to be overcome to ge nerate C 4 H 2 and hydrogen atom products. In a non gas phase system, the reactants would remain in the intermediate state due to the magnitude of the energy barrier. However in SPIAD applications, the deposition of in coming particles on the polythiophene th in films occurs in the gas phase, allowing the reactants to overcome the energy barrier. In reaction 2, the CH 2 and thiophene reactants overcome energy barriers of 3.0 kcal/mol (0.13 eV) for the site reaction and a 1.9 kcal/mol (0.083 eV) for the site reaction (Figure 7 3) The imaginary frequencies of the transition states were 1373.58 cm 1 ( site) and 1463.44 cm 1 ( site). Reactions 3 and 4 are predicted to be barrierless. In reaction 3 two highly reactive thiophene radicals generated a dithiophene molecule, C 8 H 6 S 2 Reaction 4 encompasses a thiophene and thiophene radical
112 reactants yielding the formation of a bond amongst the two molecules to produce a C 8 H 7 S 2 molecule. 7.3 Summary Computational investigations of reactions involving thioph ene and organic molecules such as CH 2 and C 2 H are carried out. Enthalpies of formations a re calculated for reactants and products. Reaction b arriers and transition states a re identified for reactions of interest. The BMK and B98 method s are shown to perfor m best in generating enthalpies of formation in comparison with experimental values. The BMK method performed best for small molecular systems while the B98 excelled in large molecular systems. Experimental results we re not available for all systems, such as thiophene radical and dithiophene molecules. Consequently, in these cases comparisons are made to the G3 and CBS QB3 methods due to their known accuracy in thermo chemistry. The B3LYP method did not execute well in comparison with the other methods. Out of the eight molecules inves tigated, the B3LYP functional i s successful in yielding enthalpy of formations close to experimental values for only two systems, C 2 H and C 2 H 2 However, this result is consistent with literature showing that B3LYP works best for small organic molecules  The resu lts from this study provide an indication of the reactivity of thiophene with organic molecules which are important for SPIAD and related polymerization and modification processes
113 Table 7 1 Comparison of the enthalpy of formation at 298K in kcal/m ol for reactant and product structures, where no charge is present. The theories used are DFT Hybrid (B3LYP, B98, and BMK), Gaussian (G2, G3), and Complete Basis Set (CBS QB3) functional methods. Compound Multiplicity B3LYP BMK B98 G2 G3 CBS QB3 Experiment C 4 H 4 S 0 33.380 27.608 30.975 29.825 27.701 23.638 27.5  C 4 H 3 S site 1 100.911 93.055 95.015 98.898 94.906 93.189 C 4 H 3 S site 1 98.095 90.353 92.226 92.610 90.736 CH 2 3 91.986 94.034 94.493 94.696 92.376 94.632 93.7  C 2 H 2 138.600 144.048 139.184 138.655 136.258 136.835 135.3 0.69  CH 3 2 32.823 35.801 36.260 35.107 33.988 35.477 35.0  C 2 H 2 1 55.847 57.107 57.623 54.888 54.014 55.471 54  C 8 H 6 S 2 1 80.309 58.523 65.369 64.193
114 Figure 7 1 O ptimized reactant structures for (from top to bottom) : Reaction 1: C 4 H 4 S + C 2 4 H 3 S + C 2 H 2 ; Reaction 2: C 4 H 4 S + CH 2 4 H 3 S + CH 3 ; Reaction 3: C 4 H 3 S + C 4 H 3 8 H 6 S 2 ; Reaction 4: C 4 H 4 S + C 4 H 3 8 H 7 S 2
115 Figure 7 2 Reaction pathway for thiophene and C 2 H forming a thiophene radical and C 2 H 2 
116 Figure 7 3 Reaction pathway for thiophene and CH 2 forming a thiophene radical and CH 3 
117 CHAPTER 8 CONCLUSIONS In order to overcome the material challenges associated with producing devices with high performance at very small length scales, it is critical to understand the fundamentals of the structure of a material and its interactions. In this dissertation, Cellular Automaton and Density Functional Theory computer simulation methods were presented to mimic a few to several behavior s of a real system, examine key aspects, and expand the understanding of organic material systems by exploring phenomena not easily accessible to experiment. Biotinylated microtubules partially coated with streptavidin gliding on surface adhered kinesin mo structures. A Cellular Automaton simulation technique is presented to model the kinetics of microtubule gliding and interactions. The analysis of the motion allows statistically meaningful data to be obtained which can be compared to experimental results, and aid in understanding the formation process of these nanostructures. With this simulation tool, multiple scenarios of interest, including those without any current experimental realization, c an be modeled quickly while varying the critical mechanistic parameters. A wealth of new insights into the assembly process has been generated. The current simulations cannot capture all aspects of the experiments. The complex shape of spools reported in t he simulations is typically not observed experimentally, because experimentally spools tend to relax to a circular state to minimize their internal stresses. Experimental observations for a direction of rotation have not be en investigated, because the chir ality of the microtubule is not modeled here.  In future work, contributions to the microtubule from kinesin motors and biotin
118 strep tavidin bonds through forces and interactions should be investigated to capture additional details concerning the formation process of nanowires and nanospools via active transport. With minor changes to the simulation tool, outcomes from microtubules bein g temporarily pinned to random surface locations can be further investigated to understand the transportation process. This simulation tool can also be applied with adaptations to other material microstructures by examining large system with a number of v ariables and transformations in recrystallization, grain growth, and phase transformations.  Manipulation in the time and space dimensions will allow the extraction of information regarding the kinetics and local variations in the microst ructure. This simulation tool can potentially also find usage in the study of biological pattern formation, such as patterns in the pigmentation of mollusc shells, [150 152] and interactions of biological cells in cell migration of systems such as actin and collagen structures, and animal herds.  Chemical reactions of thiophene with organic molecules are of interest to modify thermally deposited coatings o f conductive polymers. The reactivity of thiophene with organic molecules was examined to gain a better understanding and make predictions for surface modification applications. More specifically, the energy barriers for reactions involving thiophene and s mall hydrocarbon radicals C 2 H and CH 2 are identified. The enthalpy of formation was also calculated at 298K, and comparisons were made in the performance of DFT hybrid functional methods. The BMK and B98 functionals were found to perform best out of the se lected organic molecules in generating enthalpies of formation in comparison with experimental values. The B3LYP method did not perform well in comparison with the
119 other methods, and was only successful in yielding enthalp ies of formation close to experime ntal values for C 2 H and C 2 H 2 However, this result is consistent with literature showing that B3LYP works best for small molecules  These results from the study of the reactivity amongst thiophene with organic molecules will be important for SPIAD and related polymerization and modification processes. As a result of a lack of experimental and simulation results involving thiophene and radical systems, additional testing of more reactants is needed to examin e to develop a data base of reactions.
120 LIST OF REFERENCES  B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, P. Walter, Molecular Biology of the Cell, fourth edn., Garland Science, New York, 2002.  E. E. Glater, L. J. Megeath, R. S. Stowers, T. L. Schwarz, J. Cell Biol. 173 (2006) 545.  E. Schonteich, G. M. Wilson, J. Burden, C. R. Hopkins, K. Anderson, J. R. Goldenring, R. Prekeris, J. Cell Sci. 121 (2008) 3824.  A. Ferreira, J. Niclas, R. D. Vale, G. Banke r, K. S. Kosik, J. Cell Biol. 117 (1992) 595.  A. E. Twelvetrees, E. Y. Yuen, I. L. Arancibia Carcamo, A. F. Macaskill, P. Rostaing, M. J. Lumb, S. Humbert, A. Triller, F. Saudou, Z. Yan, J. T. Kittler, Neuron 65 (2010) 53.  P. Schwille, S. Diez, Cri t. Rev. Biochem. Mol. 44 (2009) 223.  J. Roncali, Chem. Rev. 92 (1992) 711.  R. Vidu, P. Stroeve, Ind. Eng. Chem. Res. 43 (2004) 3314.  H. Kim, J. Jang, J. Appl. Polym. Sci. 68 (1998) 785.  B. N. Kim, D. H. Lee, D. H. Han, Polym. Degrad. Stab. 93 (2008) 1214.  W. D. Hsu, S. Tepavcevic, L. Hanley, S. B. Sinnott, J. Phys. Chem. C. 111 (2007) 4199.  A. J. Heeger, J. Phys. Chem. B. 105 (2001) 8475.  Y. Choi, S. Tepavcevic, Z. Xu, L. Hanley, Chem. Mater. 16 (2004) 1924.  S. Tepavcevic, Y. Choi, L. Hanley, J. Am. Chem. Soc. 125 (2003) 2396.  S. Tepavcevic, Y. Choi, L. Hanley, Langmuir 20 (2004) 8754.  I. Rayment, H. M. Holden, M. Whittaker, C. B. Yohn, M. Lorenz, K. C. Holmes, R. A. Milligan, Science 261 (1993) 58. [17 ] A. Akhmanova, J. A. Hammer Iii, Curr. Opin. Cell. Biol. 22 (2010) 479.  H. Hess, G. D. Bachand, V. Vogel, Chem. Eur. J. 10 (2004) 2110.  R. N. Cohen, M. J. Rashkin, X. Wen, J. F. C. Szoka, Drug. Discov. Today 2 (2005) 111.
121  C. Brunner, C. Wah nes, V. Vogel, Lab Chip 7 (2007) 1263.  J. Howard, Annu. Rev. Physiol. 58 (1996) 703.  J. Howard, Mechanics of Motor Proteins and the Cytoskeleton, Sinauer Associates Publishers, Sunderland, 2001.  Y. Jeune Smith, H. Hess, Soft Matter 6 (2010) 1778.  A. Desai, T. J. Mitchison, Ann. Rev. Cell Dev. Bi. 13 (1997) 83.  R. Lahoz Beltra, S. R. Hameroff, J. E. Dayhoff, Biosystems 29 (1993) 1.  H. Hess, V. Vogel, Rev. Mol. Biotech. 82 (2001) 67.  L. Wilson, H. P. Miller, K. W. Farrell, K B. Snyder, W. C. Thompson, D. L. Purich, Biochemistry 24 (1985) 5254.  L. Song, E. J. Hennink, I. T. Young, H. J. Tanke, Biophys. J. 68 (1995) 2588.  H. Hess, C. M. Matzke, R. K. Doot, J. Clemmens, G. D. Bachand, B. C. Bunker, V. Vogel, Nano Lett 3 (2003) 1651.  A. F. Macaskill, J. E. Rinholm, A. E. Twelvetrees, I. L. Arancibia Carcamo, J. Muir, A. Fransson, P. Aspenstrom, D. Attwell, J. T. Kittler, Neuron 61 (2009) 541.  X. Wang, T. L. Schwarz, Cell 136 (2009) 163.  N. Arimura, T. Ki mura, S. Nakamuta, S. Taya, Y. Funahashi, A. Hattori, A. Shimada, C. Mnager, S. Kawabata, K. Fujii, A. Iwamatsu, R. A. Segal, M. Fukuda, K. Kaibuchi, Dev. Cell. 16 (2009) 675.  A. Kamal, G. B. Stokin, Z. Yang, C. H. Xia, L. S. B. Goldstein, Neuron 28 (2000) 449.  A. Szodorai, Y. H. Kuan, S. Hunzelmann, U. Engel, A. Sakane, T. Sasaki, Y. Takai, J. Kirsch, U. Muller, K. Beyreuther, S. Brady, G. Morfini, S. Kins, J. Neurosci. 29 (2009) 14534.  S. Niwa, Y. Tanaka, N. Hirokawa, Nat. Cell. Biol. 10 ( 2008) 1269.  C. H. Kim, J. E. Lisman, J. Neurosci. 21 (2001) 4188.  L. Guillaud, R. Wong, N. Hirokawa, Nat. Cell. Biol. 10 (2008) 19.  S. Ramachandran, K. H. Ernst, George d. Bachand, V. Vogel, H. Hess, Small 2 (2006) 330.  R. Tucker, P. Katira, H. Hess, Nano Lett. 8 (2008) 221.
122  T. Nitta, A. Tanahashi, Y. Obara, M. Hirano, M. Razumova, M. Regnier, H. Hess, Nano. Lett. 8 (2008) 2305.  H. Hess, J. Clemmens, C. Brunner, R. Doot, S. Luna, K. H. Ernst, V. Vogel, Nano Lett. 5 (2005) 629  H. Hess, Soft Matter 2 (2006) 669.  Y. Astier, H. Bayley, S. Howorka, Curr. Opin. Chem. Biol. 9 (2005) 576.  C. Brunner, H. Hess, K. H. Ernst, V. Vogel, Nanotechnology 15 (2004) S540.  T. Surrey, M. B. Elowitz, P. E. Wolf, F. Yang, F. Ne delec, K. Shokat, S. Leibler, Proc. Nat. Acad. Sci. USA 95 (1998) 4293.  T. Nitta, H. Hess, Nano Lett. 5 (2005) 1337.  P. G. Vikhorev, N. N. Vikhoreva, M. Sundberg, M. Balaz, N. Albet Torres, R. Bunk, A. Kvennefors, K. Liljesson, I. A. Nicholls, L. Nilsson, P. Omling, S. Tagerud, L. Montelius, A. Mansson, Langmuir 24 (2008) 13509.  F. Gibbons, J. F. Chauwin, M. Desposito, J. V. Jose, Biophys. J. 80 (2001) 2515.  D. V. Nicolau, D. V. Nicolau, Curr. Appl. Phys. 4 (2004) 316.  T. Nitta, A. Tanahashi, M. Hirano, H. Hess, Lab Chip 6 (2006) 881.  T. Nitta, A. Tanahashi, M. Hirano, Lab Chip 10 (2010) 1447.  D. C. Ghiglia, G. A. Mastin, L. A. Romero, J. Opt. Soc. Am. A 4 (1987) 267.  M. A. Arbib, Inform. Control 9 (1966) 177.  W. Materi, D. S. Wishart, Drug Discov. Today 12 (2007) 295.  J. P. Boon, D. Dab, R. Kapral, A. Lawniczak, Phys. Rep. 273 (1996) 55.  J. Kozicki, J. Tejchman, Granul. Matter 7 (2005) 45.  S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Numerical Mathematics and Scientific Computation, Clarendon Press, Oxford, 2001.  U. Frisch, D. D'humieres, B. Hasslacher, P. Lallemand, Y. Pomeau, J. P. Rivet, Adv. Complex Syst. 1 (1987) 649.  J. Reid, Computing 48 (1992) 219.  A. Agarwal, P. Katira, H. Hess, Nano Lett. 9 (2009) 1170.
123  S. Verbrugge, L. C. Kapitein, E. J. G. Peterman, Biophys. J. 92 (2007) 2536.  J. Gagliano, M. Walb, B. Blaker, J. Macosko, G. Holzwarth, Eur. Biophys. J. 39 (2010) 801.  G. Holzwarth, K. Bonin, D. B. Hill, Biophys. J. 82 (2002) 1784.  I. Luria, J. Crenshaw, M. Downs, A. Agarwal, S. B. Seshadri, J. Gonzales, O. Idan, J. Kamcev, P. Katira, S. Pandey, T. Nitta, S. R. Phillpot, H. Hess, Soft Matter 7 (2011) 3108.  D. G. Weiss, G. M. La ngford, D. Seitz Tutter, W. Maile, Acta Histochem. Suppl. 41 (1991) 81.  R. Kawamura, A. Kakugo, Y. Osada, J. P. Gong, Nanotechnology 21 (2010) 145603.  M. Bachand, A. M. Trent, B. C. Bunker, G. D. Bachand, J. Nanosci. Nanotechnol. 5 (2005) 718.  J. D. Crenshaw, S. R. Phillpot, H. Hess, (2011) in preparation.  D. Frenkel, B. Smit, Understanding Molecular Simulation: From Algorithms to Applications, second edn., Academic Press, San Diego, 2002.  A. V. Dobrynin, Macromolecules 38 (2005) 9304.  R. J. Stokes, D. F. Evans, Fundamentals of Interfacial Engineering, Wiley VCH, New York, 1997.  R. F. Service, Science 309 (2005) 95.  H. Q. Liu, E. D. Spoerke, M. Bachand, S. J. Koch, B. C. Bunker, G. D. Bachand, Adv. Mater. 20 (2008) 4476.  F. Pampaloni, G. Lattanzi, A. JonÂ¡Â¡, T. Surrey, E. Frey, E. L. Florin, Proc. Nat. Acad. Sci. USA 103 (2006) 10248.  M. E. Janson, M. Dogterom, Biophys. J. 87 (2004) 2723.  F. Gittes, B. Mickey, J. Nettleton, J. Howard, J. Cell Biol. 12 0 (1993) 923.  G. M. Whitesides, B. Grzybowski, Science 295 (2002) 2418.  B. Nitzsche, F. Ruhnow, S. Diez, Nat. Nano. 3 (2008) 552.  M. G. L. Van Den Heuvel, S. Bolhuis, C. Dekker, Nano Lett. 7 (2007) 3138.
124  P. J. Keller, F. Pampaloni, G. L attanzi, E. H. K. Stelzer, Biophys. J. 95 (2008) 1474.  Y. W. Gao, F. M. Lei, Biochem. Bioph. Res. Co. 387 (2009) 467.  L. Bourdieu, T. Duke, M. B. Elowitz, D. A. Winkelmann, S. Leibler, A. Libchaber, Phys. Rev. Lett. 75 (1995) 176.  J. D. Cren shaw, T. Liang, H. Hess, S. R. Phillpot, J. Comput. Theor. Nanosci. (2011).  R. Kawamura, A. Kakugo, Y. Osada, J. P. Gong, Langmuir 26 (2010) 533.  R. Kawamura, A. Kakugo, K. Shikinaka, Y. Osada, J. P. Gong, Biomacromolecules 9 (2008) 2277.  J. Clemmens, H. Hess, R. Doot, C. M. Matzke, G. D. Bachand, V. Vogel, Lab Chip 4 (2004) 83.  L. Hanley, S. B. Sinnott, Surf. Sci. 500 (2002) 500.  S. R. Kasi, H. Kang, C. S. Sass, J. W. Rabalais, Surf. Sci. Rep. 10 (1989) 1.  W. Jacob, Thin Solid Films 326 (1998) 1.  I. Luchin, Tech. Phys. Lett. 24 (1998) 159.  A. Maksimchuk, S. Gu, K. Flippo, D. Umstadter, V. Y. Bychenkov, Phys. Rev. Lett. 84 (2000) 4108.  Y. Funada, K. Awazu, K. Shimamura, H. Watanabe, M. Iwaki, Surf. Coat. Technol. 66 (1994) 514.  G. Korotcenkov, Sensor. Actuat. B Chem. 107 (2005) 209.  C. M. Chan, T. M. Ko, H. Hiraoka, Surf. Sci. Rep. 24 (1996) 1.  X. D. Xiang, X. Sun, G. Briceo, Y. Lou, K. A. Wang, H. Chang, W. G. Wallace Freedman, S. W. Chen, P. G. Schultz, Science 268 (1995) 1738.  R. S. Mane, C. D. Lokhande, Mater. Chem. Phys. 65 (2000) 1.  L. M. H. Groenewoud, G. H. M. Engbers, J. G. A. Terlingen, H. Wormeeste r, J. Feijen, Langmuir 16 (2000) 6278.  T. R. Dillingham, D. M. Cornelison, S. W. Townsend, Structural and Chemical Characterization of Vapor Deposited Polythiophene Films, third edn., AVS, Minneapolis, 1996.
125  O. Bhme, C. Ziegler, W. Gpel, Synth. Met. 67 (1994) 87.  B. Servet, G. Horowitz, S. Ries, O. Lagorsse, P. Alnot, A. Yassar, F. Deloffre, P. Srivastava, R. Hajlaoui, Chem. Mater. 6 (1994) 1809.  E. Pl, V. Hornok, D. Sebok, A. Majzik, I. Dkny, Colloid Surface B 79 (2010) 276. [102 ] F. Garnier, A. Yassar, R. Hajlaoui, G. Horowitz, F. Deloffre, B. Servet, S. Ries, P. Alnot, J. Am. Chem. Soc. 115 (1993) 8716.  G. Dennler, A. Houdayer, P. Raynaud, Y. Sgui, M. R. Wertheimer, Nucl. Instrum. Meth. B 192 (2002) 420.  J. J. Vgh, D. B. Graves, Plasma Process. Polym. 6 (2009) 320.  B. Van Der Bruggen, J. Appl. Polym. Sci. 114 (2009) 630.  E. E. Johnston, B. D. Ratner, J. Electron Spectrosc. 81 (1996) 303.  M. Ozdemir, C. U. Yurteri, H. Sadikoglu, Crit. Rev. Food Sci. 39 (1999) 457  T. J. Renk, P. P. Provencio, S. V. Prasad, A. S. Shlapakovski, A. V. Petrov, K. Yatsui, W. Jiang, H. Suematsu, P. IEEE 92 (2004) 1057.  N. J. Chou, C. A. Chang, Surface Modification of Polymers., Butterworth Heinemann Inc., Bosto n, 1994.  P. D. Prewett, Vacuum 34 (1984) 931.  M. Lorenz, V. E. Bondybey, J. Low Temp. Phys. 26 (2000) 778.  I. Jang, B. Ni, S. B. Sinnott, J. Vac. Sci. Technol. A 20 (2002) 564.  L. Hanley, Y. Choi, S. Tepavcevic, Mat. Res. Soc. Symp. Proc. 804 (2004) JJ5.3.1.  S. Tepavcevic, A. M. Zachary, A. T. Wroble, Y. Choi, L. Hanley, J. Phys. Chem. A 110 (2005) 1618.  P. Kudlacek, R. F. Rumphorst, M. C. M. Van De Sanden, J. Appl. Phys. 106 (2009) 073303.
126  M. J. Frisch, G. W. Truck s, H. B. Schlegel, G. E. Scuseria, M. A. Rob, J. R. Cheeseman, J. J. A. Montgomery, T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar, J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson, H. Nakatsuji, M. Hada M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cam mi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth, P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C. Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J. V. Ortiz, Q. Cui, A G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A. Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al Laham, C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M. W. Wong, C. Gonzal ez, J. A. Pople, Gaussian 03, Revision C.02, Gaussian, Inc., Wallingford CT, 2004.  M. Born, R. Oppenheimer, Ann. Phys. Berlin 389 (1927) 457.  C. J. Cramer, Essentials of Computational Chemistry: Theories and Models Wiley, 2004.  K. Ohno, K Esfarjani, Y. Kawazoe, Computational Materials Science: From Ab Initio to Monte Carlo Methods, first edn., Springer, Berlin, 2000.  D. S. Sholl, J. A. Steckel, Density Functional Theory: A Practical Introduction, John Wiley & Sons, Inc., Hoboken, 2 009.  S. Yeon Lee, Y. Sup Lee, Chem. Phys. Lett. 187 (1991) 302.  T. Gao, L. L. Shi, H. B. Li, S. S. Zhao, H. Li, S. L. Sun, Z. M. Su, Y. H. Lu, Phys. Chem. Chem. Phys. 11 (2009) 5124.  M. D. Wodrich, C. M. Corminboeuf, P. R. Schreiner, A. A. Fokin, P. V. R. Schleyer, Org. Lett. 9 (2007) 1851.  P. Siegbahn, J. Biol. Inorg. Chem. 11 (2006) 695.  P. E. M. Siegbahn, Q. Rev. Biophys. 36 (2003) 91.  A. D. Becke, J. Chem. Phys. 98 (1993) 5648.  A. D. Becke, J. Chem. Phys. 107 ( 1997) 8554.  A. D. Becke, Phys. Rev. A 38 (1988) 3098.  J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, C. Fiolhais, Phys. Rev. B 46 (1992) 6671.
127  A. D. Boese, J. M. L. Martin, J. Chem. Phys. 121 (2004) 3405.  J. P. Perdew, W. Yue, Phys. Rev. B 33 (1986) 8800.  J. P. Perdew, Y. Wang, Phys. Rev. B 45 (1992) 13244.  A. D. Mclean, G. S. Chandler, J. Chem. Phys. 72 (1980) 5639.  R. Krishnan, J. S. Binkley, R. Seeger, J. A. Pople, J. Chem. Phys. 72 (1980) 650.  T. Clark, J. Chandrasekhar, G. W. Spitznagel, P. V. R. Schleyer, J. Comput. Chem. 4 (1983) 294.  C. Adamo, V. Barone, Chem. Phys. Lett. 274 (1997) 242.  H. L. Schmider, A. D. Becke, J. Chem. Phys. 108 (1998) 9624.  L. A. Curtiss, K. Raghavachari, G. W. Trucks, J. A. Pople, J. Chem. Phys. 94 (1991) 7221.  L. A. Curtiss, K. Raghavachari, P. C. Redfern, V. Rassolov, J. A. Pople, J. Chem. Phys. 109 (1998) 7764.  G. P. F. Wood, L. Radom, G. A. Petersson, E. C. Barnes, M. J. Frisch, J. J. A. Montgomery, J. Chem. Phys. 125 (2006) 094106.  J. J. A. Montgomery, M. J. Frisch, J. W. Ochterski, G. A. Petersson, J. Chem. Phys. 110 (1999) 2822.  J. W. Ochterski, G. A. Petersson, J. J. A. Montgomery, J. Chem. Ph ys. 104 (1996) 2598.  G. A. Petersson, D. K. Malick, W. G. Wilson, J. W. Ochterski, J. J. A. Montgomery, M. J. Frisch, J. Chem. Phys. 109 (1998) 10570.  R. Atkinson, D. L. Baulch, R. A. Cox, J. R. F. Hampson, J. A. Kerr, M. J. Rossi, J. Troe, J. Phys. Chem. Ref. Data 28 (1999) 191.  S. G. Lias, J. F. Liebman, R. D. Levin, J. Phys. Chem. Ref. Data 13 (1984) 695.  K. Fukui, Accounts Chem. Res. 14 (1981) 363.  E. Herbst, D. E. Woon, Astrophys. J. 489 (1997) 109.  J. D. Crenshaw, S R. Phillpot, N. Iordanova, S. Sinnott, Chem. Phys. Lett. (2011) Submitted.  D. Raabe, Annu. Rev. Mater. Sci. 32 (2002) 53.
128  S. Wolfram, Nature 311 (1984) 419.  C. H. Waddington, R. J. Cowe, J. Theor. Biol. 25 (1969) 219.  D. T. Lindsay, Veliger 24 (1977) 2.  M. S. Alber, Y. Jiang, M. A. Kiskowski, Physica D 191 (2004) 343.
129 BIOGRAPHICAL SKETCH Jasmine Claren Davenport Crenshaw was born to Clarence and Segrid Davenport in Durham, NC She attended North Carolina Agric ultural and Technical State University in Greensboro, North Carolina and graduated summa cum laude with a Bachelor of Science in p hysics in May of 2005. There she was the recipient of the Nationa l Science Foundation Talent 21 s cholarship and the NC A&T Sta te University C scholarship She participated in the National Institute of Health Minority Access to Research Careers Undergraduate Student Training in Academic Research, MARC U STAR, and the Ronald E. McNair Postbaccalaureate Achievement Progra m Jasmine began her graduate study at the University of Florida in 2005 under the tutelage of Dr. Simon Phillpot In 2007, she received her Master of Science in materials science and e ngineering from the University of Florida In 201 1 Jasmine Crenshaw r eceived her Ph.D. in m aterials s cience and e ngineering from the University of Florida with a concentration in b iomaterials Her graduate education was supported by the National Science Foundation DMR 0645023 and CHE 0809376 grant s National Science Foundation Alliance for Graduate Education in the Professoriate (SEAGEP) fellowship National Consortium for Graduate Degrees for Minorities in Engineering and Science, Inc. (GEM) fellowship Eastman Kodak Company, and the University of Florida Ronald E. M cNair Graduate Assistantship