UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A microfluidic approach for antibody screening from single cells with numerical methods for modelling… Lewis, Daniel Jacob 2013

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
24-ubc_2013_fall_lewis_daniel.pdf [ 7.77MB ]
Metadata
JSON: 24-1.0073985.json
JSON-LD: 24-1.0073985-ld.json
RDF/XML (Pretty): 24-1.0073985-rdf.xml
RDF/JSON: 24-1.0073985-rdf.json
Turtle: 24-1.0073985-turtle.txt
N-Triples: 24-1.0073985-rdf-ntriples.txt
Original Record: 24-1.0073985-source.json
Full Text
24-1.0073985-fulltext.txt
Citation
24-1.0073985.ris

Full Text

A Microfluidic Approach for Antibody Screening from Single Cells with Numerical Methods for Modelling Solid-Substrate Capture Design & Validation by Daniel Jacob Lewis B.A.Sc., The University of Waterloo, 2011 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Biomedical Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) July 2013 © Daniel Jacob Lewis 2013  Abstract Antibodies are the fastest growing class of therapeutics and are also ubiquitous reagents in diagnostic and research applications. Despite the advance of many applications of antibodies in basic and applied research, the majority of novel antibody discovery and screening still typically relies on the costly, time-consuming and inefficient hybridoma approach; the development of improved approaches for the high-throughput screening and selection of monoclonal antibodies is thus of high interest. In this thesis a novel approach combining the concentration enhancement and sample manipulation benefits of polydimethylsiloxane (PDMS) multilayer microfluidics devices is validated for the high-throughput prescreening of large samples of primary antibody secreting cells. Dual functionalized beads are co-incubated with antibody secreting cells to capture both antibodies and their associated messenger RNA (mRNA) transcripts to couple sequence with affinity information in downstream fluorescence activated cell sorting (FACS). This is done in a parallel way for simultaneous cell processing over a large planar array. The device is validated with quantitative reverse transcription polymerase chain reaction (RT-qPCR) to verify successful complementary DNA (cDNA) synthesis on the bead surface and that the beads retain functional antibodies throughout the incubation and reverse transcription process. As part of this work, I also present a numerical modelling pipeline to design and characterize performance of diffusion based solid substrate cell secretion assays in microfluidic chambers prior to time-consuming fabrication and experimentation.  ii  Preface The research presented within this thesis is the result of a collaborative effort within the Hansen lab to develop high-throughput low cost methodologies for antibody screening. It builds on efforts and approaches previously developed and published by Anupam Singhal, Daniel DaCosta, V´eronique Lecault, Kevin Heyries, Oleh Petriv, and Carl Hansen. The original device designs prior to optimization were based off of those designed by Oleh Petriv, Daniel DaCosta, and V´eronique Lecault. The sieve array design was based off earlier designs by Anupam Singhal and Daniel DaCosta. The inverted PDMS fabrication protocols were adapted from work by Daniel DaCosta. The bead-based approach is the result of discussions with Carl Hansen, Kevin Heyries, and Anupam Singhal. The biotinylation of antiIgG antibodies and lysozyme was performed by Kevin Heyries. The image analysis for generating qPCR curves was performed by Mike VanInsberghe. Primer and probe sets were designed in collaboration with Mike VanInsberghe. The majority of the work done specifically to implement and test the dual-functionalized bead approach, including the design and validation, testing, experimentation, and all numerical modelling was done by the author.  iii  Table of Contents Abstract  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  Preface  List of Tables  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii 1 Introduction, Background, & Motivation . . . . . . . . . . . 1.1  Antibodies  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4  1.1.1  Structure and Diversity . . . . . . . . . . . . . . . . .  4  1.1.2  Hybridoma Method for Antibody Screening and Selection  1.2  1  . . . . . . . . . . . . . . . . . . . . . . . . . .  6  Microfluidics . . . . . . . . . . . . . . . . . . . . . . . . . . .  9  1.2.1  Multilayer Soft Lithography (MSL)  9  1.2.2  Polymerase Chain Reaction (PCR) on Microfluidic Devices  1.2.3  . . . . . . . . . .  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  12  Research Statement . . . . . . . . . . . . . . . . . . . . . . .  14  2 Device Design . . . . . . . . . . . . . . . . . . . . . . . . . . . .  16  1.3  2.1  Microfluidics for Antibody Screening  11  Prototype Design 2.1.1  . . . . . . . . . . . . . . . . . . . . . . . .  19  Device Fabrication . . . . . . . . . . . . . . . . . . . .  25  iv  Table of Contents 2.2  MATLAB-COMSOL Pipeline . . . . . . . . . . . . . . . . . .  25  2.2.1  . . . . . . . . . . . . . . . . . .  28  3 Results & Discussion . . . . . . . . . . . . . . . . . . . . . . . .  30  Model Configuration  3.1  Modelling of Capture Efficiency  . . . . . . . . . . . . . . . .  30  3.2  Prototype Design Optimization  . . . . . . . . . . . . . . . .  41  . . . . . . . . . . . . . . . . . . . . . .  43  3.2.1 3.3  3.4  3.5  Button Valves  Methods & Assay Development  . . . . . . . . . . . . . . . .  50  3.3.1  Bead Functionalization  . . . . . . . . . . . . . . . . .  50  3.3.2  Protocols . . . . . . . . . . . . . . . . . . . . . . . . .  52  Bead Testing . . . . . . . . . . . . . . . . . . . . . . . . . . .  64  3.4.1  Validation in Tube Reactions with qPCR . . . . . . .  64  3.4.2  Antibody Binding & Imaging . . . . . . . . . . . . . .  69  Validation on Microfluidic Devices with qPCR  . . . . . . . .  73  4 Conclusions & Recommendations . . . . . . . . . . . . . . . .  80  Works Cited . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  83  Appendices A Numerical Modelling Pipeline . . . . . . . . . . . . . . . . . . A.1 Introduction  . . . . . . . . . . . . . . . . . . . . . . . . . . .  A.1.1 Finite Element Method and COMSOL Background  96 97  . 100  A.1.2 Calculating Flux - 1D Heat Diffusion Example . . . . 104 A.1.3 Lagrange Multipliers for the Finite Element Method . 108 A.2 Mean & Total Intensity Modelling . . . . . . . . . . . . . . . 115 B Derivations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 B.1 Concentration in a Box . . . . . . . . . . . . . . . . . . . . . 120 B.2 Analytical Model for Perfect Absorber Bead Capture  . . . . 127  C Materials & Methods  . . . . . . . . . . . . . . . . . . . . . . . 131  C.1 LabView Interfaces  . . . . . . . . . . . . . . . . . . . . . . . 131  v  Table of Contents D Coding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 D.1 Bead and Cell Chamber Configuration  . . . . . . . . . . . . 133  D.2 Bead Capture Model Construction . . . . . . . . . . . . . . . 141 D.3 Bead Capture Data Analysis . . . . . . . . . . . . . . . . . . 151  vi  List of Tables 3.1  Reverse transcription (RT) reaction mixture used both for in tube reactions (at half volume) and on-device (as shown) . .  3.2  PCR reaction mixture with TaqMan  ®  probe, used both for  in-tube reactions (at half volume) and on-device (as shown) . 3.3  58 63  Polymerase chain reaction (PCR) reaction mixture with EvaGreen® , used both for in-tube reactions (at half volume) and on-device (as shown) . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.4  63  Conditions and cycle threshold values for the various validation reactions carried out in tubes to support the validity of the dual-functionalized bead approach . . . . . . . . . . . . .  66  3.5  Ct values for Figure 3.24 . . . . . . . . . . . . . . . . . . . . .  67  3.6  Antibody binding assay with varying protein surface coverage and either direct or indirect capture of specific antibodies . .  3.7  70  Assay for bead binding ability with 100% surface coverage of lysozyme against anti-HEL with and without a heat lysis step prior to reverse transcription . . . . . . . . . . . . . . . . . .  3.8  Testing detection limit for functionalized beads against labelled antibodies . . . . . . . . . . . . . . . . . . . . . . . . .  3.9  70 71  Results for the antibody binding reaction conditions shown in Figure 3.8 . . . . . . . . . . . . . . . . . . . . . . . . . . .  71  A.1 Comparison of basis function derivatives vs. Lagrange multiplier values to compute boundary flux (with different numbers of finite elements for the 1D heat example) . . . . . . . . . . 110  vii  List of Figures 1.1  Flowchart describing the antibody screening methodology . .  3  1.2  The five subclasses of immunoglobulin . . . . . . . . . . . . .  5  1.3  IgG antibody structure . . . . . . . . . . . . . . . . . . . . . .  6  1.4  Hybridoma technology method . . . . . . . . . . . . . . . . .  8  1.5  Illustration of the multilayer soft lithography fabrication procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  1.6  Illustration showing the structure and actuation of multilayer soft lithography pneumatic valves . . . . . . . . . . . . . . . .  2.1  17  Illustration depicting the surface bi-functionalization of the streptavidin coated polystyrene microbeads . . . . . . . . . .  2.3  11  Flowchart illustrating the required steps a functional device prototype must enable for successful methodology validation .  2.2  10  18  Basic unit cell of the device design with parallel stochastic loading inverted array illustrating the functional steps required for capture, co-incubation, lysis, and reverse transcription . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  22  2.4  Illustration showing the operation of the button valve design  23  2.5  AutoCAD layout file showing the inlet tree and valving structure enabling individual access to four subarrays for comparative reaction conditions and loading schemes within the same device . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  2.6  24  Flow chart illustrating the numerical modelling pipeline for bead capture of released molecules in a microfluidic chamber  27  viii  List of Figures 2.7  Graphical representation of an example for each of the four main bead distributions used for model construction as generated by the MATLAB configuration script . . . . . . . . . .  3.1  Surface flux of 140 individual beads in a gradient configuration within a cuboid chamber of side-length 150 µm . . . . . .  3.2  29  31  Percentage of released mRNA molecules from a lysed cell captured on adjacent beads vs. incubation time in a microfluidic chamber as modelled by the numerical modelling pipeline . .  3.3  Percentage of released molecules captured by the beads over time for the various bead configurations . . . . . . . . . . . .  3.4  34  Scatter plot of total captured molecules on each individual bead in a chamber vs. the linear distance from the cell . . . .  3.5  32  37  Bar graph indicating the total number of beads included in a mean or total intensity calculation vs. linear distance from the cell for 60 beads in each of the four configurations . . . .  3.6  38  Bar graph indicating the total number of beads included in a mean or total intensity calculation vs. linear distance from the cell for 100 beads in each of the four configurations . . . .  3.7  39  Bar graph indicating the total number of beads included in a mean or total intensity calculation vs. linear distance from the cell for 140 beads in each of the four configurations . . . .  3.8  40  Brightfield image showing chamber contents moving towards one side in the direction of flow, with rectangular capture chambers . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  3.9  41  Isosurface plot of fluid flow velocity for a rectangular chamber design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  3.10 Streamlines for rectangular box inverted chamber geometry .  43  3.11 Example bright-field image from the volume containment button valve testing . . . . . . . . . . . . . . . . . . . . . . . . .  44  3.12 Results of button valve permutation testing for chamber volume retention . . . . . . . . . . . . . . . . . . . . . . . . . . .  45  3.13 Flow results of Figure 3.12 with SPR lanes removed . . . . .  46 ix  List of Figures 3.14 Rank order results from Figure 3.13 . . . . . . . . . . . . . .  47  3.15 Isosurface plot for a circular chamber design . . . . . . . . . .  48  3.16 Streamlines for cylindrical inverted chamber geometry . . . .  49  3.17 Flowchart illustrating the validation protocol . . . . . . . . .  52  3.18 Example cell loading distribution . . . . . . . . . . . . . . . .  56  3.19 Composite image of the nuclear and cytoplasmic stain for imaging lysis buffer effects on live cells (BaF3 cells shown) . .  57  3.20 One dimensional diffusion model . . . . . . . . . . . . . . . .  59  3.21 Graphical depiction of the solution to the concentration profile in a one dimensional box of constant concentration reagent inflow as shown in Equation 3.4 . . . . . . . . . . . . . . . . .  60  3.22 Image of ROX fluorescent dye loaded device with insufficient flow time to achieve uniform reagent concentration across the array . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  61  3.23 qPCR results for in-tube bead validation . . . . . . . . . . . .  65  3.24 Comparison of different reverse transcription temperatures and heat lysis with bead’s surface functionalized with 50% biotinylated goat anti-mouse IgG and 50% b2m reverse primer 67 3.25 Effect of proportional protein surface coverage on the RTPCR efficiency of in-tube reactions at equivalent concentrations to device . . . . . . . . . . . . . . . . . . . . . . . . . .  68  3.26 Example 350 ms exposure of concentrated beads originally incubated at 1 bead per nL with labelled antibody at 4 nM concentration showing the ease of distinguishing signal from background . . . . . . . . . . . . . . . . . . . . . . . . . . . .  72  3.27 Individual chamber qPCR curves from on-device PCR with imaging at each annealing step . . . . . . . . . . . . . . . . .  74  3.28 Example of manually curated heat-map for cell-count occupancy of each chamber in the array . . . . . . . . . . . . . . .  75  3.29 qPCR cycle threshold value heat map for the device with occupancy shown in Figure 3.28 . . . . . . . . . . . . . . . . .  76  3.30 Superimposed heat-map of cell occupancy count . . . . . . .  77  x  List of Figures 3.31 Box plot of cycle threshold value vs. occupancy number of each chamber in the array . . . . . . . . . . . . . . . . . . . .  78  3.32 Agarose gel electrophoresis of RT-PCR products . . . . . . .  79  A.1 Diagram for the analytical single bead model of a perfect absorber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  97  A.2 Concentration in a well-mixed solution over time with a single 5µm diameter bead acting as a perfect absorber in a 2nL chamber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  98  A.3 Arbitrary geometry in two dimensions for the coefficient form PDE with domain Ω in the plane R2 . . . . . . . . . . . . . . 100 A.4 Number of elements vs. flux accuracy for the traditional finite element method . . . . . . . . . . . . . . . . . . . . . . . . . . 111 A.5 Depiction of the geometry used for the basic diffusion and capture model. Drawing units are in µm . . . . . . . . . . . . 112 A.6 Cell and bead flux vs. time for the model depicted in Figure A.5 for the derivative-based approach to calculating boundary flux . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 A.7 Cell and bead flux vs. time for the model depicted in Figure A.5 for the Lagrange multiplier approach to calculating boundary flux . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 A.8 Percent difference between the theoretically predicted measured mean and total intensity between a low-end number of beads (60) and a high-end number of beads (140) for different configurations in a 150 µm square chamber . . . . . . . . . . . 116 A.9 Mean intensity measurement both theoretical and predicted actual measurement (based on a threshold of detection of 5000 molecules) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 A.10 Total intensity measurement both theoretical and predicted actual measurement (based on a threshold of detection of 5000 molecules) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 B.1 System Definition . . . . . . . . . . . . . . . . . . . . . . . . . 120  xi  List of Figures B.2 Diagram for the analytical single bead model of a perfect absorber . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 B.3 Solution concentration vs. time for analytical bead capture model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 C.1 LabView interface used for timing and solenoid/valve control. 131 C.2 Main LabView interface used for microscope control and imaging. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132  xii  Acknowledgements First and foremost I want to acknowledge and thank my supervisor, Dr. Carl Hansen, for letting me play around with microfluidics in his lab and for teaching me about the field and the intersection of biology and engineering. I also want to thank the thesis committee members Dr. Steve Withers and Dr. Karen Cheung for their critical reading and assessment. I wish to thank everyone whom I worked with in the Hansen lab for making my two years at UBC an educational and enjoyable experience; your hard work and dedication have been inspiring during my time here. I would particularly like to thank Kevin, Mike, and Adam for their helpful discussions and advice throughout the course of this project. Kevin has been incredibly patient and helpful with regards to interpreting experimental results and designing the next one. Mike was very generous to donating his time to the image analysis of some difficult to handle images, and for helping design the GAPDH and b2m primer/probe assays for validation work. I also am thankful to my family and friends for their support and encouragement and for making me go out now and again to relax.  xiii  Chapter 1  Introduction, Background, & Motivation The development and optimization of a microfluidics based pre-screening platform for identification and selection of targeted antibodies is investigated, with ultimate applications in novel antibody discovery for therapeutics, diagnostics, and research. Current antibody screening techniques largely rely upon the dated hybridoma method which is costly, time consuming and inefficient. Microfluidic approaches offer a novel way to screen for targeted antibodies without these shortcomings. A microfluidic platform and assay protocol are designed that may be easily scaled to allow for highthroughput, rapid, and parallel screening and selection from tens to hundreds of thousands of primary B lymphocytes (B cells), reducing the cost and time required to identify targeted antibodies of interest. The system is validated on murine hybridoma cells (fused mouse B cells and myelomas). In addition, numerical methods and software systems are developed to reduce the design time and number of iterative steps required in the microfluidic design process for solid substrate assays. The main thrust of this thesis is to use microfluidic compartmentalization to implement an antibody display technology in which both functional antibodies and the genetic sequence that codes for them are linked together on microparticles. An overview of the designed screening methodology is shown in Figure 1.1. We aim to ultimately clone antibody genes from highaffinity targeted antibodies into expression vectors for recombinant production in mammalian cell lines. Both affinity and sequence information is thus required from each single antibody secreting cell. To this end, single  1  Chapter 1. Introduction, Background, & Motivation cells and dual functionalized micro-beads are co-incubated in microfluidic chambers allowing for quick capture of secreted antibodies and their associated messenger RNA (mRNA). Reverse transcription (RT) is performed to link complementary DNA (cDNA) onto the bead surface. Subsequently, the beads are eluted from the device and co-incubated with fluorescently labelled antigen. Beads containing antibody specific to the labelled antigen are then detected and isolated by fluorescence activated cell sorting (FACS) into single wells of PCR plates. From there, polymerase chain reaction (PCR) can be used to amplify the cDNA off the bead surface, followed by DNA sequencing to obtain the transcript sequence information. From there, the antibody genes would be cloned into expression vectors for recombinant production in mammalian cell lines. The work herein describes a design validation of the pre-screen assay used to couple the affinity and sequence information onto the surface of the micro-beads.  2  Chapter 1. Introduction, Background, & Motivation  Figure 1.1: Flowchart describing the antibody screening methodology. Single cells and dual functionalized microbeads are co-incubated in microfluidic chambers allowing for quick capture of secreted antibodies and their associated messenger RNA (mRNA). Reverse transcription (RT) is performed to link complementary DNA (cDNA) onto the bead surface. Subsequently, the beads are eluted from the device and co-incubated with fluorescently labelled antigen. Beads containing antibody specific to the labelled antigen will respond fluorescently to excitation light and are sorted using fluorescence activated cell sorting (FACS) into single micro-wells. From there, polymerase chain reaction (PCR) can be used to amplify the cDNA off the bead surface, followed by DNA sequencing to obtain the transcript sequence information. 3  1.1. Antibodies  1.1  Antibodies  Antibodies are globular proteins of the vertebrate immune system (referred to as immunoglobulins or Ig), produced by B cells to defend against foreign invaders by recognizing and binding to foreign antigen [1, 2]. Their functions vary by subclass and include: neutralization and blocking; enabling effective macrophage opsonization, phagocytosis, and degranulation; and activation of the classic complement pathway [1]. Antibodies can have extremely high affinity to their target moiety, ranging from 10−7 -10−11 M [1]. The high affinity and specificity of antibodies makes them extremely useful many biological and biomedical applications including: affinity capture, protein characterization, and diagnostics [3–9]. Their specificity, excellent tolerability, and low toxicity also make antibodies the fastest growing class of therapeutic, with a global market of over $50 billion and a yearly growth rate in excess of 10% [3]. Despite their widespread utility, traditional approaches to antibody screening and selection remain expensive, inefficient, and time-consuming [9–11].  1.1.1  Structure and Diversity  Immunoglobulins are sub-divided into five major subclasses based on structure and biological effects (Figure 1.2): IgD, IgE, IgG, IgA, and IgM [1]. Of these, IgG is the most often used in biological applications [1, 2]. Each class of antibody has an identical constant region (Figure 1.3) that is conserved within species [12]; this can be used as a common capture site for monoclonal or polyclonal anti-IgG antibodies and is used as such in this assay. The structure of IgG forms a symmetric Y-shaped molecule, comprised of 4 separate chains (Figure 1.3). At the outer end of each protein ‘arm’ is a variable region (identical on both sides of the protein) that serves as the antigen binding site. These 4 chains are made up of two identical ‘heavy’ and two identical ‘light’ chains (Figure 1.3). The heavy chains form a constant Fc region at the base of the body. Both the heavy and light chains form a variable Fv region in the arms each contributing 3 complementarity determine regions (CDR) to each arm of the antibody — for a total of 4  1.1. Antibodies  Figure 1.2: The five subclasses of immunoglobulin (Reproduced with permission from Martin Br¨ andli). six CDRs. Both these chains must be sequenced for successful recombinant expression. In humans, the heavy chain sequence is created by a random recombination of variable (V), diversity (D) and junction (J) regions during B cell maturation, producing ∼9000 unique combinations [1, 13]. Similarly, the light chain sequence is produced by a recombination of Vκ regions and Jκ regions or Vλ and Jλ regions, producing ∼200 unique combinations [2, 14]. Following combinatorial selection of heavy and light chains, this yields nearly 2 million unique combinations. Diversity is further increased by trimming and addition of bases [14] giving rise to ∼ 1014 possible unique molecules. While useful, this diversity makes primer design difficult.  5  1.1. Antibodies  (a) Basic IgG antibody structure showing the symmetrical Y shape and variable region antigen binding sites on the ‘arms’. (Public domain image).  (b) Heavy and light chain IgG antibody structure, showing constant and variable regions. Two disulfide bonds connect the heavy chains and another two connect the heavy chain to the light chain on each arm for a total of 4 disulfide bonds per molecule. (Public domain image).  Figure 1.3: IgG antibody structure.  1.1.2  Hybridoma Method for Antibody Screening and Selection  The broad adoption and use of antibodies in research and targeted therapeutic applications to date has been enabled by the development of hybridoma technologies [15]. Most antibodies selected for human use are selected by the hybridoma method from engineered mice or rats expressing human IgG. However, individual B cells do not produce antibody at a sufficient rate to allow for pre-screening of specificity or affinity by conventional techniques [14]; antibody producing B cells also have a short life-span, only surviving for several days [1]. Therefore, to obtain a measureable source of antibody, B cells are fused with myelomas (a cancerous B cell) that are selected for their ability to generate viable fusions with B cells. The hybrid cell is referred to as a hybridoma [16]. Notably, the availability of a suitable fusion partner is only established for rodents - mice, rats, and rabbits. This restricts antibody selection by hybridoma approaches to these species. A large number 6  1.1. Antibodies of B cells must be hybridized and each cell allowed to proliferate in culture for over a week to obtain sufficient quantities of antibodies in culture supernatant to allow for screening1 [17, 18]. Figure 1.4 illustrates the hybridoma process from rodent immunization to multiplication and harvesting. While the vast majority of antibodies produced for therapeutic and research use have been derived from some variant of the hybridoma method, the method has several shortcomings [15, 18, 19]. Extracted B cells from immunized rodents will typically number around 1,000,000 cells; however, the myeloma fusion process successfully generates only ∼1000 stable clonal lines per extraction of which 100 to 200 may be antigen specific depending on the antigen [14, 15, 18, 19]. This implies that only a small portion of antibody diversity is ultimately tested, resutling in a loss of many high-value antibodies. By enabling the direct analysis of antibodies secreted from single cells, microfluidic approaches offer a novel way to screen for targeted antibodies without these shortcomings.  1  Screening is typically done with immunoprecipitation or enzyme-linked immunosorbent assays (ELISA)  7  1.1. Antibodies  Figure 1.4: Hybridoma technology method. In (1), a rodent is immunized against an antigen of interest. In (2) B cells are harvested from the spleen of the animal and myeloma cells are extracted and cultivated in (3). The B cells and myelomas are fused in (4) to produce hybridoma lines. In (5) these hybridoma lines are isolated into single cells and grown in culture. Step (6) is screening of the antibodies produced from these monoclonal lines against the antigen of interest. Those producing relevant antibody are kept (6a) and the rest are discarded (6b). In step (7) the cells are expanded in vitro (7a) or in vivo (7b). Finally, the antibodies are harvested in (8). (Reproduced with permission from Martin Br¨andli).  8  1.2. Microfluidics  1.2  Microfluidics  The term microfluidics refers to technologies used to manipulate fluids and their suspended components at the micro-scale, typically involving channel dimensions in the tens of microns and fluid volumes in the nanolitre to femtolitre range [20]. Starting from the first micro-machined silicon-based devices focused on chromatography and ink-jet printing [21–23], microfluidics has emerged as a broad and distinct field incorporating a wide range of fabrication techniques and materials including silicon, glass, photoresist and elastomers — with a wide range of applications in chemistry and biology [20, 24–26]. Most prominently, devices fabricated with the elastomer polydimethylsiloxane (PDMS) using soft lithography approaches [27, 28] have emerged as the dominant medium for microfluidic applications in biology. This is due to the several advantages of PDMS such as its mechanical properties (enabling fast, simple and re-usable valves), transparency to visible light, biocompatibility and gas permeability, along with the ability to rapidly and cheaply iterate through prototypes when compared to silicon or glass-based fabrication protocols [20, 29–31]. The combination of small volumes and complex fluidic control have made the handling of individual cells and small quantities of molecules in a highthroughput parallelized fashion much easier and more precise than traditional approaches [20, 30, 32–36]. In addition, the low reagent consumption and concentration enhancement have drastically reduced experimental cost and time for many applications when compared to typical bench-top protocols [20, 37].  1.2.1  Multilayer Soft Lithography (MSL)  The breakthrough of microfluidics technologies from early proof-of-concept to the forefront of biology has been enabled by the rapid prototyping technique called multilayer soft lithography and the invention of the stacked layer valve [27, 28]. Multilayer soft lithography is essentially a replica molding technique using the elastomer PDMS. Used as a two part viscous fluid, 9  1.2. Microfluidics a base and a curing agent are mixed and poured over a micro-fabricated master mold whose surface structure is the inverse of the desired channel geometry of the microfluidic device. The mixed fluid cures at room temperature in a process called room-temperature vulcanization (RTV). From there, the hardened PDMS device can be removed and the master mold reused. A graphical illustration of the multilayer soft lithography technique is shown in Figure 1.5.  Figure 1.5: Illustration of the multilayer soft lithography fabrication procedure. A blank silicon wafer is spin-coated with photoresist. The photoresist is exposed to ultraviolet light, cross-linking negative photoresist. Areas unexposed are subsequently etched by submersion in developer solvent. This process is repeated to build up the inverse three dimensional structure of the microfluidic device. Following fabrication and parylene coating (to prevent irreversible PDMS adhesion to the surface), a PDMS mixture of base and cross-linking agent is poured over the surface and cured via RTV. Multiple layers may be bonded via oxygen plasma treatment of the surface to build up the final structure [38–40]. Access ports are punched to allow for external fluidic connection and the entire structure is bonded via oxygen plasma treatment to a glass slide. 10  1.2. Microfluidics The second major innovation allowing for complex microfluidic devices was the invention of computer controlled pneumatic valves, integrated into the device structure [28]. A thin control layer below a thick flow layer, or a thick control layer above a thin flow layer forms an intersecting channel structure separated by a thin membrane. When pressurized, the control layer membrane deflects into the flow layer, restricting flow; this is shown in Figure 1.6. The low Young’s modulus of PDMS allows the membrane to be deflected with actuation pressure of 20-30 psi, and may be repeated indefinitely.  Figure 1.6: Illustration showing the structure and actuation of multilayer soft lithography pneumatic valves. The lower control layer is pressurized to deflect a thin membrane up into the thick flow layer, restricting flow. Note that the red thin PDMS membrane is actually part of the control layer; it has been given a separate colour for clarity.  1.2.2  Polymerase Chain Reaction (PCR) on Microfluidic Devices  Real-time quantitative polymerase chain reaction (qPCR), and more specificity, reverse transcription qPCR (RT-qPCR) is a method for reverse transcribing and amplifying the cDNA of RNA transcripts while quantitatively monitoring the increase in molecule quantity by means of probes [41], molec-  11  1.2. Microfluidics ular beacons [42], or intercalating dye [43]. It relies upon PCR to amplify reverse transcribed DNA from as little as a single molecular copy [44]. An RNA transcriptase first reverse transcribes sample RNA into a cDNA chain which is subsequently amplified exponentially by a DNA polymerase through multiple rounds of denaturation and extension using thermocycling. Differences in starting quantity are measured by monitoring the reaction and determining how many cycles of PCR are required to pass through a predetermined threshold of detection. This value, referred to as the cycle threshold value (Ct ), is used to quantify the amount of initial cDNA present relative to other samples, and thus the number of transcripts initially in solution. The total number of molecules N present in solution after n cycles of PCR is related to the initial number of copies ci and an amplification efficiency factor A by: N = A · ci · 2n  (1.1)  This quantity of molecules N will be linearly related to the measured fluorescence by a variety of factors including the dye used, the length of the DNA molecule, and its sequence. Identical samples with different treatments can be compared relative to each other in order to assess total workflow efficiency. Or, by comparison with samples of known concentration and similar amplification efficiency, the initial copy number can be calculated in an absolute sense. In this thesis, RT-qPCR is used to assay the success of reverse transcription and linking of heavy and light chain (or endogenous controls in proof of concept experiments) cDNA to the dual functionalized bead surface.  1.2.3  Microfluidics for Antibody Screening  Recently a number of different approaches for single-cell isolation and screening via cloning (bypassing the hybridoma fusion process) have been developed [45–47]; single antibody secreting B cells are isolated from tissue samples with FACS, using labelled antibodies against B cell surface markers. The antibody heavy and light chains are then amplified using single-cell reverse transcription polymerase chain reaction (RT-PCR) followed by cloning 12  1.2. Microfluidics into expression vectors for recombinant expression in bacterial, yeast, or mammalian cells [14, 45–47]. Following culture and recovery, the antibodies are subsequently screened for binding affinity and specificity to the desired target antigen. This process has been successful but still requires laborious cloning and cell culture prior to screening. Microfluidics is an enabling technology well-adapted to single-cell manipulation and analysis. When working with finite amounts of analyte, as is the case with single cells, the small volume confinement of microfluidics provides a large concentration enhancement that makes it possible to directly analyze the antibodies generated from a single B cell. The small volumes, and therefore concentration enhancement, allows for the design of devices and protocols that can isolate and screen individual B cells without the need for hybridoma fusion approaches or the laborious cloning and culture required for screening in other approaches. By comparison, the incubation of single cells in a traditional micro-well plate (typically ≥10 µl in volume) would require several weeks of incubation before reaching concentrations detectable by ELISA [14]. The microfluidic chambers used in this thesis have a volume of approximately 1 nl, allowing for the accumulation of nanomolar concentrations of antibody from a single B cell in approximately 20 min. Several approaches focused on exploiting this small volume effect using micro-fabricated systems have been developed using micro-wells [48–50], emulsions [51], and a combination of the two [52]. The micro-engraving approach by Love et. al. [48] uses a planar array of micro-wells in a PDMS substrate with manual pipetting of single antibody secreting cells into each well. Following gravity settling, the array is sealed against a glass slide functionalized with anti-IgG antibodies to capture secreted antibodies. The array is then separated from the glass slide and the slide treated with fluorescently labelled antigen to identify spots with relevant antibodies. The cell corresponding to that spot can be manually removed from the microwell array for subsequent expansion. Jin et. al [49] use a similar micro-well array to isolate individual antibody producing cells; however, the surface of the array is coated with anti-IgG antibody. As antibody from the antibody secreting cells diffuse up out of the well it is captured radially on the surface 13  1.3. Research Statement outside the well. Introduction of fluorescently labelled antigen causes a ring of fluorescence around wells producing antibody with affinity to the target antigen. These cells can then be manually extracted. K¨oster et. al. [51] designed a droplet-generating microfluidic device for single-cell encapsulation of hybridoma cells for discovery of antibodies that inhibit enzymatic activity. Finally, DeKosky et. al. [52] have developed an approach capturing single antibody secreting cells in micro-wells followed by emulsion linkage PCR to obtain the heavy and light chain sequences. This would then be followed by cloning and recombinant expression then screening. While these approaches have shown the ability to screen through up to tens of thousands of single cells, they are still largely at the proof-ofconcept stage and each has drawbacks which limit their functional utility. The droplet-based approach makes reagent addition and removal difficult and complicates the recovery of cells and sequences following screening. The linkage PCR method requires cloning prior to screening, and the microwell approaches, while successful, have only shown the ability to identify a handful of antigen-specific antibodies out of those screened and have labour intensive manual pipetting steps. The work described herein is focused on the development and validation of an alternative high-throughput approach that may easily be scaled to prescreen hundreds of thousands of single antibody secreting cells in a rapid and parallel fashion, reducing the cost and time required to identify targeted antibodies of interest. Following selection of cells and sequencing of the heavy and light chains, the genes may be cloned and expressed for further kinetic screening.  1.3  Research Statement  Despite their utility, methods for the low-cost and high-throughput screening of antibody secreting cells for antigen-specific antibodies remain elusive. The specific goals for the research outlined in this thesis are: 1. Design, develop, and validate a microfluidic device and assay for pre14  1.3. Research Statement screening of targeted antibodies from single antibody secreting cells 2. The device should be easily scalable, robust, and amenable to automation 3. Design and code a numerical modelling pipeline to enable the characterization and optimization of microfluidic designs for solid substrate capture in micro-wells  15  Chapter 2  Device Design The microfluidic device designed must enable the screening methodology described in Figure 1.1. The specific requirements are laid out in sequential order in Figure 2.1. The device should enable high-throughput and rapid capture and isolation of single cells in microfluidic chambers followed by washing of the isolated cells to remove any cellular debris and free RNA from lysed cells as well as any antibody in the supernatant. After the cells are washed the device must enable loading of beads (dual functionalized with anti-IgG antibody and reverse primer) into the same chambers without losing the cells or causing lysis. The chambers must then be individually isolated to allow for incubation and capture of secreted antibody on the bead surface. Following incubation, the device must allow for cell lysis without any cross-talk between capture chambers. After cell lysis, the mRNA must be captured on the same beads, followed by introduction of RT mix followed by synthesis of cDNA off the primer on the bead surface. Finally, the device must not only allow for elution and recovery of the dual functionalized beads for downstream sorting but also allow for on-device qPCR for RT and cDNA validation. The primary utility of the bead-based approach lies in the inherent coupled information it provides: the antibody specificity information through the introduction of labelled antigen and subsequent FACS; and the associated sequence information for that antibody through the bound cDNA following reverse transcription. This is achieved with a dual surface functionalization of both anti-IgG antibodies targeting the Fc region for the species being studied, and sequence specific or multiplexed primers2 . Multi2  As an alternative approach, beads functionalized directly with the antigen of interest were also investigated; following specific binding to the bead, an Fc specific labelled  16  Chapter 2. Device Design  Figure 2.1: Flowchart illustrating the required steps a functional device prototype must enable for successful methodology validation. plexed reactions can be performed on a single bead, and multiple transcripts copied simultaneously as each bead has enough binding sites to theoretically capture the entire transcriptome of a single cell [53]. Polystyrene beads surface coated with streptavidin were used to couple biotinylated primers and anti-IgG antibodies to the surface for mRNA capture for cDNA synthesis and capture of secreted IgG (See Figure 2.2). To ensure functional utility of a validated design beyond the prototype stage, there are two additional key requirements to consider: first, the device must easily scale to enable sample processing at the throughput required to sample a large amount of the antibody repertoire (∼100,000 cells); second, the device design should be amenable to automated integration in routine operation. In addition to device architecture design and optimization for loading and flow control, a numerical modelling pipeline was also designed to exanti-IgG antibody would be introduced prior to FACS. Both approaches were seen to be effective.  17  Chapter 2. Device Design  Figure 2.2: Illustration depicting the surface bi-functionalization of the streptavidin coated polystyrene microbeads. The large blue sphere represents the polystyrene bead, the red spheres represent surface-linked streptavidin molecules, the green spheres represent biotin — affixed to the primers and anti-IgG antibodies. Note that the diagram is not drawn to scale. amine the capture of secreted antibodies and release of mRNA following lysis and how bead number, distribution and chamber geometry affected the capture. It is difficult to accurately control these variables experimentally and to achieve consistent precise results for resulting fluorescence (in the case of antibody capture) and impossible to measure mRNA capture directly. There are two direct cases of interest for which a well-developed bead capture model is useful: first, how long a lysed cell must be incubated with a given number of beads in a chamber to capture the majority of mRNA molecules of interest released during lysis; second, what is the surface coverage of antibody on a given bead in the chamber over time, from a nearby secreting cell, and how does this vary with distribution and number of beads? To answer these questions a MATLAB-COMSOL finite element method based numerical modelling pipeline was developed. The pipeline  18  2.1. Prototype Design was designed such that others could use it with little to no experience in programming or modelling software.  2.1  Prototype Design  The isolation of single cells in a high-throughput manner requires a parallel loading architecture. This requirement in combination with the needed ability to wash and flush-in additional reagents led to a novel device design incorporating an inverted planar array, based off previously published work for microfluidic cell culture arrays [54]. An inverted chamber design allows for parallel stochastic loading of both cells and subsequently beads into the device and — after settling — enables washing and additional reagent to be brought into the chamber without removing the contents. Co-incubation of the cells and beads in these chambers allows for antibody capture and — following lysis — of released mRNA. The ability to flush reagents allows for reverse transcription mix to be brought into the chambers for synthesis of cDNA on the bead surface. Finally, the device may simply be inverted and flushed to release the beads for subsequent off-device screening. A diagram of a device unit cell illustrating the aforementioned steps is shown in Figure 2.3 . The bead-based approach enables coupling of both the heavy and light chain information with the antibody specificity information allowing for simultaneous elution of all samples prior to FACS (See Section 3.4.1, Section 3.3.2 and Section 1). However, while the design enables the collection of specificity information, it is essentially a binary assay, and downstream kinetics measurements would be essential to proper antibody selection. The primary advantage of a stochastically loaded parallel architecture like this over a serial loading design is the increase in throughput it enables. With near-instant priming of the entire device, loading of the cells and beads across many thousands of chambers can be achieved in several minutes. The initial design strategy for isolated lysis was a heat lysis protocol. The device chambers could be locked down with the isolation valves and the entire device placed on a thermocycler at elevated temperature (87 °C for 7 19  2.1. Prototype Design min) to permeabilize the cells and release the mRNA. This allows for a simple device architecture and a guarantee of no cross-talk between the chambers. However, it was discovered that heat-lysis would not be a satisfactory lysis method due to reduced RT efficiency and antibody denaturation (See Section 3.4.1). This motivated the need for a mild chemical lysis technique. Flowing a chemical lysis buffer over the array as designed would result in cells near the inlet lysis first and spreading their contents throughout the rest of the array, before the entire device could be primed. To avoid this, a novel approach incorporating push-down “button” valves over cylindrical chambers was designed . This design allows for the array to be primed with a lysis solution while keeping individual chambers isolated. Following priming, the cells can be individually exposed to the lysis solution to release the mRNA without cross-talk. The cylindrical chamber was given a lip of greater or equal diameter connecting it to the lateral flow channels. Above the chamber a circular valve was created which seals off the inner chamber by line-contact with the edge when actuated, while still allowing flow through the surrounding lip; thus the array can be flushed with additional reagents without any inflow or outflow from the chambers themselves, preventing cross-contamination or early lysis. An illustration of this design is shown in Figure 2.4. The upper disk size and chamber depth were selected such that upon opening the button valve, a solution loaded at 10x concentration would be diluted to 1x upon mixing with the chamber volume. Prototype devices designed for validating the bead approach are small in scale (8 devices per 4 inch silicon wafer). However, as per the design requirements, the end-goal is a large device with tens to hundreds of thousands or more of individual chambers. The design of the device is such that it is easily expanded in size to make large devices while also lending itself to easy automation. The bottom of the device — below the cell capture chambers — is a thin PDMS membrane bonded directly to glass, allowing for high quality images useful for image analysis in bead and cell counting. The pneumatic valves are all controlled through custom LabView software and could be run from a prescheduled program; multiple small devices could 20  2.1. Prototype Design also theoretically be run in parallel to achieve the throughput desired. The prototype device was designed with a tree-based valving structure for inlet lines and outlet lines that allows isolated access to four individual subarrays within the large array with input from one of four separate reagent lines. This allows for multiple reactions and loading schemes to be compared simultaneously on the same device (See Figure 2.5).  21  2.1. Prototype Design  Figure 2.3: Basic unit cell of the device design with parallel stochastic loading inverted array illustrating the functional steps required for capture, coincubation, lysis, and reverse transcription. Cells are stochastically loaded (1) and allowed to settle (2). The cells are washed and uncaptured cells are removed from the array lines (3). Beads are then stochasticaly loaded (4), allowed to settle (5), and washed (6). The large valve in the center is a novel button valve design allowing for chamber isolation while allowing the array to be flushed. The secondary isolation valves above the flow lines are actuated to isolate each chamber allowing for incubation and exposure to reagent in the flow lines (after releasing the button valves) without crosstalk. In (7) the button valve is actuated allowing the array to be primed with lysis solution. Actuating the secondary isolation valves and opening the button valve allows for diffusive mixing of the lysis buffer with the chamber contents, bursting the cell (8).  22  2.1. Prototype Design  Figure 2.4: Illustration showing the operation of the button valve design. When not pressurized, the device functions as a simple inverted array with cylindrical chambers, with most of the flow passing over the deep capture well allowing for cell and bead capture (6). When actuated, the button valve isolates the capture chamber with line contact around the edge of the chamber, preventing inflow or outflow (7) and allowing the array to be primed with additional reagents. The presence of the circular lip allows for flow around the chamber during actuation. After this priming, the adjacent secondary isolation valves can be actuated for isolation and the button valve released to expose the chamber contents to the new reagents (8).  23  2.1. Prototype Design  Figure 2.5: AutoCAD layout file showing the inlet tree and valving structure enabling individual access to four subarrays for comparative reaction conditions and loading schemes within the same device. Dark blue indicates SPR photoresist fabricated channels for restricted flow during valve actuation. Light blue indicates SU-8 photoresist fabricated channels for the hydration lines and upper chamber lip. Green indicates the cell-capture chambers. Red indicates the valving layer.  24  2.2. MATLAB-COMSOL Pipeline  2.1.1  Device Fabrication  Devices were fabricated primarily using standard photolithography and multilayer soft lithography techniques [27, 28]. Masks for each different height were created by laser etching photoresist covered, chrome coated quartz mask blanks. These masks were used to progressively build up the inverse mold of the device using SPR and SU-8 photoresist with a separate mold for the control valve layer. Once the mold was complete it was parylene coated using vapour deposition to create a 1-100 nm layer of parylene, to prevent PDMS adhesion. From there, a 10:1 ratio of PDMS and cross-linking agent was mixed, degassed, and poured onto the control layer to form the thick component of the device. The flow layer containing the inlet lines and inverted chambers was fabricated by spin-coating a thin layer of PDMS onto the mold followed by an inverse PDMS stamping technique. The two layers were then bonded together via oxygen plasma bonding and access holes were punched to allow for external tubing connection [38–40]. The device is finally affixed to a glass slide again by oxygen plasma bonding.  2.2  MATLAB-COMSOL Pipeline  COMSOL is a finite element method (FEM) multi-physics modelling language and environment [55]. It is used to: define geometries and their associated boundary and initial conditions, automate finite element mesh construction, automate basis function construction, and automate derivation of the stiffness, mass, and loading matrices used for solving the basis function coefficients. Details on the mechanics of COMSOL can be found in the various user manuals [56–60]. Further explanation of the FEM than is described herein can be found in [61–66]. MATLAB® is a language and interactive environment for numerical computation [67] and is integrated into COMSOL through the LiveLink™ interface. The link allows for the construction and manipulation of models through MATLAB and then sends the resulting configuration file to COMSOL for matrix construction and computation.  25  2.2. MATLAB-COMSOL Pipeline The pipeline designed uses MATLAB as the front-end to collect user input for the geometry and modelling parameters; this input is then used to create parameter and coordinate text files for the chamber geometry and bead and cell placement. These text files are pulled into a second script to create the model and feed it to the COMSOL kernel for processing. The user can create as many coordinate and parameter files as desired and the pipeline will process each sequentially and independently without additional user input. The output from each constructed model is saved in text files which are then opened by MATLAB for automatic data processing and visualization. Data saved in these files includes the flux over each bead at each instant in time modelled by the discretization along with its position in the chamber and position relative to the cell, all the parameters used for the construction of the model, and the concentration everywhere in the chamber at each instant in time so the user is able to construct additional data analysis approaches if desired. The flowchart in Figure 2.6 illustrates the process. Full coding details for the chamber configuration can be found in Appendix D.1, for the model construction in Appendix D.2, and for the data analysis in D.3. Full details on the modelling mathematics and design can be found in Appendix A.  26  2.2. MATLAB-COMSOL Pipeline  27  Figure 2.6: Flow chart illustrating the numerical modelling pipeline for bead capture of released molecules in a microfluidic chamber.  2.2. MATLAB-COMSOL Pipeline  2.2.1  Model Configuration  The user is first prompted to input basic geometric parameters of the design to be modelled including chamber dimensions and shape, number of beads (and distribution if desired), and the size and number of cells (and position if desired). This information is used to first construct an area over which the beads and cells may be positioned; from there, the user specified number of components and distribution types are used to randomly generate locations for each bead and cell in the chamber subject to the imposed constraints. An example visual output for four different distribution types for a chamber with one cell against a wall and 100 beads is shown in Figure 2.7. The code for this step is contained within Appendix D.1. From there, the process proceeds as illustrated in Figure 2.6. The four example configurations shown were chosen based on patterns observed experimentally for the distribution of beads within a chamber. Early designs used rectangular capture chambers and the lack of a fluid vortex at the bottom of the chamber (See Figure 3.10 and Figure 3.16) meant that when reagent swapping occurred, the chamber occupants would often be pushed towards one side of the chamber in the direction of flow (Represented by the ‘gradient’ and ‘close’ configurations). Other work in the lab has investigated the possibility of separation of the cells and beads onto opposite sides of the chamber as shown in the ‘far’ configuration, and for comparison, a completely random configuration was also included in the analysis. For each configuration type and number of beads, 3 random distributions were created and sent through the modelling pipeline. The results shown in Section 3.1 represent an average of the data produced by those models.  28  (b) ‘Far’ configuration  (c) ‘Gradient’ configuration  (d) ‘Random’ configuration  29  Figure 2.7: Graphical representation of an example for each of the four main bead distributions used for model construction as generated by the MATLAB configuration script. There are 100 beads of 5µm diameter, the cell is 10µm in diameter, and the chamber floor is a 150µm by 150µm square.  2.2. MATLAB-COMSOL Pipeline  (a) ‘Close’ configuration  Chapter 3  Results & Discussion 3.1  Modelling of Capture Efficiency  As described in Section 2.2.1 and Appendix A.1.3, the modelling pipeline employed is able to accurately keep track of flux over each element of a Dirichlet boundary within the domain being discretized and thus track the number of molecules that would be present on the surface of each bead. To illustrate this, an example output graphic from a gradient distribution containing 140 beads is shown in Figure 3.1. The colormap spans red to blue, with red indicating a higher surface flux, and blue a lower surface flux. The first question posed is “how long must a lysed cell be incubated with a given number of beads in a chamber to capture the majority of mRNA molecules released after lysis?” This was modelled with a cylindrical capture chamber of 100 µm radius, 125 µm deep, and with only 15 beads present in the chamber, randomly distributed alongside a single cell containing 1000 mRNA transcripts. The model assumed that these molecules would be released uniformly over a 20 s period during which the permeabilization of the plasma membrane occurred. A conservatively chosen diffusion constant of 4 × 10−11 m2 s−1 was used, calculated from the radius of gyration of a 1200 base-pair long RNA molecule [68, 69]. The results of this model are shown in Figure 3.2. As is seen in Figure 3.2, nearly 100% of all released molecules are captured within 15 min of lysis. To ensure that nearly all molecules are captured, when accounting for variability in diffusion, lysis, and bead distribution, a 2x safety factor was used to choose the 30 min incubation time used when developing the protocol discussed in 3.3.2. When using this 30 min incubation time in practise, little to no cross-contamination was observed 30  3.1. Modelling of Capture Efficiency  Figure 3.1: Surface flux of 140 individual beads in a gradient configuration within a cuboid chamber of side-length 150 µm. Beads are shown coloured with blue representing low surface flux and red representing high surface flux. The cell is in wireframe. and nearly every chamber with intact cell occupation resulted in qPCR ‘hits’ indicating the model accuracy was sufficient to enable design and protocol decisions. The second question posed is “how the surface coverage of antibody on a given bead in a chamber varies with distribution and number of beads?” A related question that follows from this is more to the point: what is the optimal number of beads in a chamber, and their optimal distribution, to achieve accurate and quantitative imaging results and proper FACS sorting? The pipeline was set-up with the configuration replicates as described in  31  3.1. Modelling of Capture Efficiency  Figure 3.2: Percentage of released mRNA molecules from a lysed cell captured on adjacent beads vs. incubation time in a microfluidic chamber as modelled by the numerical modelling pipeline. Section 2.2.1. The flux over each bead in the chamber in each configuration at each instant in time was saved to file and used to construct the results. Several outstanding questions in the lab related to this antibody capture problem were: whether mean or total intensity of the beads in a chamber represented a more accurate and quantitative approximation of the production rate of antibodies from an antibody secreting cell (See Appendix A.2); how long you must incubate for adequate bead capture; and what is the capture efficiency of a given configuration. The main factor contributing to differences in these cases, assuming that each experiment produces a similar distribution of beads in each of the chambers would be the number of beads in a chamber; more beads would effectively distribute the flux over a larger area and thus the mean intensity would be expected to be lower for an otherwise identical cell; however, with more beads, there is a greater likelihood of 32  3.1. Modelling of Capture Efficiency a molecule being captured at any given moment in time, and thus one might expect this to compensate for the lowered average due to quantity. But how does the intensity of an arbitrary bead in the chamber vary with time? For different bead configurations, the results can be surprisingly different. The capture efficiency is easily compared between configurations by looking at the total percent capture of molecules released from a secreting cell over time. Figure 3.3 shows that all numbers of beads and configurations rapidly approach high capture efficiency, with the ‘far’ configuration lagging behind the others. It was noted that a sufficient quantity of antibody on the average bead required for imaging was achieved after only 20 min of incubation. The main point of note is that total capture efficiency is relatively similar irrespective of the number of beads in the chamber, so a closer look at the differences in individual bead variability is required.  33  3.1. Modelling of Capture Efficiency  34 Figure 3.3: Percentage of released molecules captured by the beads over time for the various bead configurations. Results shown for 60 beads, 100 beads, and 140 beads. The dashed line represents the theoretical maximum where the entire floor of the chamber is a perfectly absorbing surface.  3.1. Modelling of Capture Efficiency To compare individual bead intensity, a scatter plot of the total captured molecules on each individual bead as a function of their linear distance from the cell, at a given time point is plotted (2 h of incubation time). This is shown in Figure 3.4; notably, for all but the far configuration, it is primarily the beads closest to the cell which take up most of the flux, at the expense of the other beads. With the far configuration, the concentration has time to build up across the chamber such that when the molecules reach the beads on the other side of the chamber, many more of the beads see an identical flux. This causes more beads to be included throughout the measurement space and the intensity is comparable from one bead to the next. There are three main considerations regarding this distributed flux: for in-situ comparisons, more beads included in the calculation averages out technical noise such as chip defects and lens variability; for sorting via FACS, it is desired for every bead that came from the same chamber be as close in intensity as possible to accurately sort beads from chambers relative to other chambers; and finally, it should be noted that the beads are not perfectly homogeneous. By maintaining a greater number of beads at an elevated intensity, the negative effect of bead heterogeneity is greatly reduced as there is likely to be at least one high-performing bead within the chamber that has captured enough antibody to be properly representative of the antibody specificity. Greater homogeneity in signal intensity, assuming relatively homogenous bead surfaces and sizes, is created with fewer numbers of beads in the chamber, as illustrated by Figure 3.5, Figure 3.6 and Figure 3.7; with sufficient surface coverage necessary for imaging achieved in all configurations. Furthermore, this is clearly advantageous for coupling sequence and antibody specificity information as fluorescent beads without sequence information, or vice versa result in a loss of information and time. Fewer beads are thus highly preferable for capture. In addition, it is noted that for the ‘far’ configuration, the effect is maintained even for higher number of beads (Though requiring longer incubation times, as shown in Figure 3.3). Future work should investigate the ability to more accurately control the number and distribution of beads on the chamber floor to take full advantage of these 35  3.1. Modelling of Capture Efficiency effects. For the purposes of this thesis, beads were used in quantities of 5-20 beads per chamber, randomly distributed following stochastic loading.  36  3.1. Modelling of Capture Efficiency  37 Figure 3.4: Scatter plot of total captured molecules on each individual bead in a chamber vs. the linear distance from the cell.  3.1. Modelling of Capture Efficiency  38 Figure 3.5: Bar graph indicating the total number of beads included in a mean or total intensity calculation vs. linear distance from the cell for 60 beads in each of the four configurations.  3.1. Modelling of Capture Efficiency  39 Figure 3.6: Bar graph indicating the total number of beads included in a mean or total intensity calculation vs. linear distance from the cell for 100 beads in each of the four configurations.  3.1. Modelling of Capture Efficiency  40 Figure 3.7: Bar graph indicating the total number of beads included in a mean or total intensity calculation vs. linear distance from the cell for 140 beads in each of the four configurations.  3.2. Prototype Design Optimization  3.2  Prototype Design Optimization  The device was initially designed with a chamber depth of 100 µm, a width of 100 µm, and a length (along the flow channel direction) of 150 µm. The valves were designed with a width of 150 µm, and a membrane thickness of 15 µm to actuate over flow channels 20 µm in height. During preliminary experimental investigation into device performance for capture and washing, it was noted that even with low inlet pressures of 1-3 psi, that the contents of the inverted chambers would shift along the bottom of the chamber in the direction of flow during reagent swapping and flushing steps, resulting in a crowding of chamber contents up against the far wall of the chamber (See Figure 3.8); this made manual bead and cell counting difficult, and would likely make future image analysis based automation in cell-counting unreliable. At higher flow rates (∼5-9 psi) — required for reagent swapping in a reasonable time-frame — some of the chamber contents would be pulled back up into the flow stream and either exit the device, or re-deposit in downstream chambers.  Figure 3.8: Brightfield image showing chamber contents moving towards one side in the direction of flow, with rectangular capture chambers. To investigate and approach a design solution to this problem, the flow characteristics of the capture chambers were modelled in COMSOL, a multiphysics modelling software suite. A simulation of the device chamber revealed that even at low flow rates, a significant degree of flow occurs down into the chamber. The chamber was re-designed from a depth of 100 µm to 41  3.2. Prototype Design Optimization a depth of 150 µm to reduce the flow penetration at the bottom of the chamber. Figure 3.9 shows that at this depth, most of the flow passes through the upper half of the chamber.  Figure 3.9: Isosurface plot of fluid flow velocity for a rectangular chamber design. Chamber dimensions are 150 µm deep, 100 µm wide, 100 µm in length with an upper flow channel 20 µm in depth. The chamber was modelled with an inlet pressure of 5 psi, no back pressure, and a no-slip boundary condition on the walls. The increased depth was effective at preventing loss of chamber contents at high flow rates, however, the chamber contents were still pushed to one side. Figure 3.10 shows that, despite the majority of flow occurring through the top-half of the chamber, streamlines are still present along the bottom floor of the chamber, causing a shift of contents towards the wall. This was observed to occur until the beads and cells were immediately adjacent to the wall, where no flow occurs.  42  3.2. Prototype Design Optimization  Figure 3.10: Streamlines for rectangular box inverted chamber geometry. Modelled with an inlet pressure of 5 psi, no-slip boundary conditions on the walls, and no back pressure. Streamlines penetrate to the bottom of the chamber causing a lateral shift in chamber contents, even at low pressure.  3.2.1  Button Valves  A permutation array of various ratios using both SPR and SU-8 photoresist for rounded and flat channels respectively was created to test different designs. The flow channel and upper chamber lip were fabricated at 15 µm in depth with the flow channel at 100 µm in width. The button valves were actuated at 25 psi with an inlet pressure on the flow-stream of 7 psi. The membrane thickness for valving was 15 µm. The array was initially filled with food dye coloured water, followed by actuation of the button valves. The array was then continually flushed with deionized water to remove the food dye. The intensity of each chamber was monitored over a period of 3 h to assess inflow of water into the dye filled chambers. The chambers whose intensity remained constant and thus retained their volume, were successful designs. The results of this testing are shown in Figure 3.12 and an example bright-field image is shown in Figure 3.11. Flow testing (Figure 3.12) indicated that SPR formed chambers com43  3.2. Prototype Design Optimization  Figure 3.11: Example bright-field image from the volume containment button valve testing (Figure 3.12). pletely collapsed with valve actuation allowing no flow throughout the length of the flow line and thus no volume could be replaced. The SU-8 lanes however were effective for certain inner diameter to ring ratios. Removing the SPR results and re-plotting the table shows a clear pattern (Figure 3.13). Plotting the results in rank order of increasing outer ring radius to inner chamber diameter ratio produces the chart shown in Figure 3.14. It is clear that for ratios greater than 1.7, volume retention is successful.  44  3.2. Prototype Design Optimization  45  Figure 3.12: Results of button valve permutation testing for chamber volume retention. Colours are indicated by legend in the lower left. Flow direction was along columns. The Top row indicates whether any flow was possible in the column below. The material indicates the type of photoresist to fabricate the upper lip of the chamber and thus whether it was rounded or flat, respectively. The H/D ratio on the right side indicates the chamber height to diameter ratio. All chambers on this array were 80 µm in height. Valving membrane thickness was 15 µm. The numbers on the left indicate the diameter of the chamber. The numbers along the bottom indicate the outer ring (lip) radius. The numbers within the array indicate the outer ring radius to inner chamber radius ratio.  3.2. Prototype Design Optimization  46  Figure 3.13: Flow results of Figure 3.12 with SPR lanes removed.  3.2. Prototype Design Optimization  Figure 3.14: Rank order results from Figure 3.13 Once a suitable button valve design was identified, modelling and testing was done to optimize the chamber width relative to the flow line as well as chamber depth. Experimental results indicated that a capture chamber with diameter equal to or greater than the flow channel width was most effective at capturing cells as the majority of the flow passed directly over the capture chamber. Any smaller meant that many cells that would otherwise have been captured were stuck on the lip outside of the chamber and were subsequently washed away. COMSOL modelling of flow profiles indicated similar penetration depth into the chamber as the rectangular chamber profile (See 47  3.2. Prototype Design Optimization Figure 3.15).  Figure 3.15: Isosurface plot for a circular chamber design. Chamber dimensions are 125 µm deep, 100 µm in diameter, with an upper disk diameter of 170 µm, and a flow channel 20 µm in depth. Chamber was modelled with an inlet pressure of 5 psi, no back pressure, and a no-slip boundary condition on the walls. Streamline investigations of various cylindrical designs showed that for chambers whose diameter was equal to that of the upper flow channel width, a slow-moving circular vortex of flow would develop along the bottom of the chamber (See Figure 3.16). The vortex — which was not present using the rectangular cross-section chambers — creates a flow profile resembling a fluid cushion that prevents excessive flow penetration along the bottom of the chamber. If the modelling was accurately representative of the flow characteristics of a real device, it would be expected that even at elevated inlet pressures, the lower chamber contents should remain at the bottom of the chamber. A cylindrical chamber array was fabricated using the design constraints identified in Figure 3.12 and it was observed that, even at inlet pressures of up to 7-15 psi that the chamber contacts would remain within the chamber. This validated the modelling approach and simultaneously solved the lateral shift problem that made cell counting difficult. Moreover, 48  3.2. Prototype Design Optimization it was observed that cells and beads remained essentially stationary in these chambers, showing that the magnitude of the flow within the vortices is also significantly attenuated relative to the rectangular chambers.  Figure 3.16: Streamlines for cylindrical inverted chamber geometry. A circular fluid vortex is created in the bottom of the chamber, providing a slow-moving fluid cushion to prevent streamline penetration. The chamber was modelled with an inlet pressure of 5 psi, no back pressure, and a no-slip boundary condition on the walls. A new device was thus designed with circular chambers and button valves where the inner chamber and flow channel are 100 µm in diameter, and the outer lip is 170 µm in diameter. The height and spacing was chosen such that a lysis buffer loaded at 10x concentration, would dilute to 1x concentration upon opening the button valve and exposing the two volumes. Adjacent to the button valves (See Figure 2.4) are standard valves for restricting flow. When the button valves are opened to expose cells to the lysis buffer, the adjacent secondary isolation valves are closed to prevent cross contamination via diffusion from the lysed cells. The volume of the array can thus be replaced without removing the contents of the inverted chambers. Initial device designs incorporated large chambers similar to those used  49  3.3. Methods & Assay Development for inverted cell culture arrays. While necessary for effective cell culture, these are not needed for efficient reverse transcription and the large unitcell spacing creates problems with chip dehydration; the gas permeability of PDMS becomes a disadvantage with water vapour escaping at elevated temperatures, changing the salt and other reagent concentration of reactions, or drying out the chambers completely; these losses can range from minimal to over 90% fluid loss. Unexpectedly, this dehydration during the RT incubation step resulted in a strong vacuum being formed, where the chambers were visibly deformed inwards. This resulted in the button valves being permanently “stuck” closed and prevented recovery of beads from the device. To avoid this problem, the array was re-designed with smaller, more closely spaced chambers for higher density and more cells per device. This reduces dehydration during long incubation and reverse transcription steps by maintaining a higher effective water density in the core of the device, slowing the evaporative process. The new design reduced dehydration in the device core, and should theoretically be sufficient based on past observations. In practise however, dehydration was still observed spreading radially inward from the outermost chambers of the array where the effective water density is lower. To combat this, the valve channels over the array — used to isolate the chambers and normally filled with Krytox oil — were filled with reverse transcription solution to maintain a source of osmotically neutral solution and increase the local water density in an attempt to offset outward diffusion of water. Hydration lines were also added to each side of the planar array. This, in combination with continually flowing solution after the introduction of button valves made array flushing possible, reduced dehydration to a minimum.  3.3 3.3.1  Methods & Assay Development Bead Functionalization  Polystyrene microbeads of 5 µm diameter and surface-linked with streptavidin (Bangs Laboratories) were coupled to the antigen, anti-IgG antibody,  50  3.3. Methods & Assay Development and primers via biotin-streptavidin bonds. Biotin-Streptavidin affinity is one of the strongest non-covalent molecular interactions known. It is stable at elevated temperatures and essentially irreversible. Thus, following binding, any surface modifications are also essentially irreversible, allowing for downstream processing without losing the coupled information these beads provide. The beads used from Bangs Laboratories have ∼2,000,000 binding sites per bead. However, due to the large size of antibodies, steric hindrance will prevent all sites from being occupied. In addition, the close proximity of primer to the surface and to adjacent antibodies may make it difficult for the reverse transcriptase and DNA polymerase enzymes to begin synthesis. For this reason various molar coverages were tested, and a spacer molecule was used to separate the primers from the bead surface. The reverse strand primer for the transcript of interest is covalently biotinylated at the 5’-end (Integrated DNA Technologies) to allow for the 5’-3’ synthesis activity of reverse transcriptase along with a triethylene glycol (TEG) spacer to reduce steric hindrance near the bead surface; similarly, antibodies are ordered which have been covalently biotinylated (abcam). The antibodies and primers are mixed together at the desired molar ratio and then co-incubated for several hours with the streptavidin-coated microbeads (at 10x the molar amount of surface sites available to ensure 100% surface coverage, in an attempt to avoid non-specific adsorption to the surface), followed by several washing steps. The end result is a functionalized microbead with tuneable relative surface coverage of protein and primer. This is illustrated in Figure 2.2. After functionalization, the beads are re-suspended to a stock concentration of 30 beads/nl in pH 8.0 Tris-EDTA (TE) buffer spiked with 0.1% Tween-20 and 0.05% BSA. It was found that the addition of both these reagents was necessary to prevent bead self-aggregation and sticking to the tube walls, even when Lo-Bind® tubes (Eppendorf) were used. In addition, the BSA was observed to improve the efficiency of the RT reaction, possibly due to the molecular crowding effect.  51  3.3. Methods & Assay Development  3.3.2  Protocols  An iterative process involving experimental investigation and modelling was used to design a robust device. This section is written so as to guide the reader through the experimentation protocol for a working device and include the information used to arrive at the protocol as it is written. This work focused on achieving consistent on-device RT-qPCR of chosen transcripts from bi-functionalized beads. The general approach used for the on-device work is illustrated in Figure 3.17.  Figure 3.17: Flowchart illustrating the validation protocol.  Device Priming with Solution The entire device is first ‘primed’ with Dulbecco’s phosphate buffered saline (PBS, Gibco) containing 1% bovine serum albumin (BSA) prior to loading beads or cells (12 psi) for several minutes. The BSA coats the channel walls and significantly reduces cell and bead adhesion when compared to control loading experiments (Confirmed by visual inspection). Typically, PDMS-based microfluidic devices are also primed with a 0.1% Tween-20 solution which has shown to be effective at reducing RNA and DNA adhesion to channel walls, promoting efficient RT and PCR down to the single molecule level [34, 70]. However, it was noticed that when Tween20 was used for priming that false positive ‘hits’ would be observed in the qPCR curves for empty chambers suggesting partial cell permeabilization was occurring prior to lysis. For this reason Tween-20 was removed from all solutions ‘seen’ by the cells until cell lysis and bead capture of released mRNA was complete. 52  3.3. Methods & Assay Development Bead Loading Following channel priming, the device is flushed with PBS (12 psi) for several minutes to remove any unbound BSA. Unwashed lines would exhibit BSA aggregation/crystallization causing irreversible clogging. Beads which have been resuspended3 to a concentration of 10-30 beads/nl in PBS + 0.05% BSA are then prepared for loading by rapid pipetting in the tube to prepare a well mixed solution free of aggregates. Beads are loaded into the desired sub-arrays (6-8 psi). The ability to selectively load beads and subsequently wash the lines without affecting other sub-arrays allows for the use of multiple bead types in different sub-arrays over the course of a single experiment. After the respective sub-arrays are primed with beads, flow is stopped by closing the outlet valve. The beads are then allowed to settle to the bottom of the chambers. Beads are loaded prior to cells as it was found that — due to the lack of Tween-20 present in the priming solution — some mRNA from permeabilized cells within the cell suspension would adhere to the walls; the beads would pick up this contamination as they were loaded and empty chambers would show qPCR ‘hits’. This effect is not observed if the beads are loaded prior to cell loading. It is important that cells and beads are allowed to fully settle into chambers prior to washing steps to avoid losing chamber contents. A quick calculation using Stoke’s Law for terminal velocity (Equation 3.1 [71]) with water density of 1000 kg m−3 (ρf ) [72], bead density of 1050 kg m−3 (ρp ) (Bangs manufacturer data), gravity of 9.8 m s−2 (g), dynamic viscosity of 0.001 N s m−2 (µ) [72], a bead radius of 2.5 µm (R) and a distance to travel of 150 µm to the bottom of the chamber gives a settling time of 3.68 min. This was verified to be accurate experimentally, and in practise the beads are given 5 min to settle to the bottom of the chamber. After settling, the subarrays are flushed with PBS to remove beads residing within the flow Re-suspension is performed by centrifuging at 10 000 rpm (11 500 m s−2 reactive centrifugal acceleration) for 3 min followed by removal of the supernatant 3  53  3.3. Methods & Assay Development channel and inlet lines. vs =  2 (ρp − ρf ) 2 gR 9 µ  (3.1)  Cell Loading For validation experiments, HyHel-5 cells were used. HyHel-5 is a murine hybridoma cell line that expresses IgG specific to the model protein hen egg lysozyme (HEL). After bead loading and flushing, HyHel-5 cells were removed from suspension culture and washed several times in cold PBS + 0.05% BSA4 . The cells are re-suspended after the final wash to a concentration of 1 × 106 cells ml−1 in cold PBS + 0.05% BSA. Cells are then loaded onto the chip at 1-3 psi. Once the array has been primed the outlet port is shut and the cells above chambers allowed to settle to the bottom for 5 min. The cell loading process is stochastic and related to the loading concentration and array geometry. If it assumed that the cells are well-mixed and loading events are independent, then cell loading density will follow a Poisson distribution (Equation 3.2 [73]). An example loading distribution is shown in Figure 3.18. Despite attempts to maintain a constant concentration of cells across experiments there was significant variability in observed chamber occupancy; achieving consistent cell loading without trapping with cell traps [34] or sequential valve isolation with serial loading [35] remains an open problem but was not addressed in the current work. Collating cell occupancy data from several experiments yields Poisson distributions with λ varying from 0.03 to 15 . The lower number gives, for a 1000 chamber array, approximately 29 chambers with single cells, and likely 0 for the rest of the array. For higher loading densities, this gives 367 chambers with single cells, but 183 chambers with two cells, 61 chambers with three cells, 15 chambers Washed by centrifuging at 1500 rpm for 5 min (1725 m s−2 reactive centrifugal acceleration) removing supernatant, and re-suspended by gentle pipetting 5 In (3.2) f is the frequency of event X = k where k = 0, 1, 2, ... and λ is the expected value of X for a given distribution. 4  54  3.3. Methods & Assay Development with four cells, and 3 chambers with five cells. f (k; λ) = P r(X = k) =  λk e−λ k!  (3.2)  55  3.3. Methods & Assay Development Figure 3.18: Example cell loading distribution. The numbers indicate the quantity of cells present in each chamber. Only sub-array 1 and sub-array 4 were primed with cells. 56  3.3. Methods & Assay Development Incubation and Lysis After settling, the arrays are again flushed with PBS (containing 1.5 mM Mg2+ ) to remove cells residing in the flow channels; this is done until no visible cells or cellular debris remain in the device outside the chambers (12 min). PBS is used for washing to maintain an osmotic balance between the chamber volume and interior of the cell, preventing lysis or permeabilization prior to isolation. Following flushing, the button valves are actuated at 25 psi to create line contact between the upper membrane and the edge of the cylindrical chamber, isolating the co-incubating cells and beads from the flow channel while still allowing the array to be flushed with reagent. A 10x lysis buffer consisting of 10% Triton X-100 and 20 mM MgSO4 is loaded onto the array. A device was designed to test the effect of various lysis buffers on cell membrane integrity. Cells were stained with a nuclear stain and a cytoplasm stain to allow for imaging of lysis (Figure 3.19). Various lysis buffers were tried; it was found that Triton X-100 at a 1% concentration was effective at permeabilizing the plasma membrane of mammalian cells without disrupting the nucleus. This is important since the release of genomic DNA can inhibit RT and PCR reactionson device in small volumes and/or lead to sticking of cells and beads to the device surface.  Figure 3.19: Composite image of the nuclear and cytoplasmic stain for imaging lysis buffer effects on live cells (BaF3 cells shown). After loading, the secondary isolation valves are actuated and the button valves de-actuated to allow for diffusive mixing of the lysis buffer into the cell chamber; this diffusive mixing is observed to happen within 2 min of 57  3.3. Methods & Assay Development de-actuation and the cell is visibly permeabilized. Following lysis, the cell is incubated with the beads for 30 min at 55 °C by securing the device to a thermal cycler with a flat plate. Modelling pulsed release by a cell (contact with lysis buffer) with adjacent beads shows greater than 95% capture after a 15 min incubation at elevated temperature (See Section 3.1). The elevated temperature was also found to improve mRNA binding efficiency, which is likely due to disruption of secondary structure. Reverse Transcription After incubation, the chip temperature is reduced to 10 °C prior to loading RT reaction mix. The RT reaction mix is shown in Table 3.1. Reagent RNase-free Water 5x First Strand Buffer dNTP Mix SuperScript™ III RT Tween 1% BSA 1% RNase Out Total Volume  Volume (µl) 22.5 10 2.5 2.5 5 5 2.5 50  Table 3.1: Reverse transcription (RT) reaction mixture used both for in tube reactions (at half volume) and on-device (as shown). The 10 °C chip temperature and minimum 1 mM divalent magnesium concentration present at all times following cell lysis are chosen to drastically limit the ability of the RNA-primer complex to dissociate during reagent flushing of the array [74, 75]. Prior to flushing the chambers with RT mix, the button valves are again closed and the array flushed with RT mix; this is done to remove the viscous 10% Triton X-100 solution from the array. Once the device is primed with RT mix, the button and secondary isolation valves are de-actuated and the full device — including chambers — is flushed with RT mix (6-9 psi). While the mRNA lost during reagent swapping should be negligible, it 58  3.3. Methods & Assay Development is still desired to minimize the time required for chambers to be exposed to one another to further reduce any chances of cross-contamination between chambers in addition to shortening the required experiment duration. Imaging experiments with food-dye-loaded water indicate the upper flow channel rapidly fills with incoming solution which then diffuses down into the deep cell capture chambers. While some molecular concentration is lost along the length of the array due to sequential diffusion into each chamber, the chambers also fill faster than a purely diffusive model would predict due to convection down into the chamber itself. Therefore, a diffusion model assuming constant concentration at the top of the chamber can be used to approximate the required filling time for a chamber of given dimensions. A simple one dimensional model is constructed as shown in Figure 3.20. The governing equation for a diffusive system in one dimension is given by Fick’s second law of diffusion (3.3): ∂φ ∂2φ =D 2 ∂t ∂x  (3.3)  x=0  x=L  φ(0, t) = C0  ∂φ =0 ∂x  Figure 3.20: One dimensional diffusion model for constant concentration diffusing into a chamber. A simplified model for predicting flow-time required for reagent swapping across the linear array. Dimensions and boundary conditions are shown with initial condition of φ(x, 0) = 0. The solution to the system as depicted in Figure 3.20 and described by Equation 3.3 is shown in Equation 3.4 (For full derivation details see Appendix B.1). A graphical representation of this solution is shown in Figure 3.21 for a diffusion constant of 3 × 10−11 m2 s−1 (chosen conservatively to  59  3.3. Methods & Assay Development represent a slow moving large protein such as the reverse transcriptase in solution). To achieve greater than 90% replacement of volume within the chambers, this model predicts a flow time of at least 15 min is required — assuming a sufficiently high flow rate to maintain a high concentration within the flow channel. ∞  φ(x, t) = C0 + n=1  4(Ci − C0 ) −(2n − 1)2 π 2 Dt (2n − 1)πx exp sin π(2n − 1) 2L 4L2 (3.4)  In (3.4), φ(x, t) is the concentration of species φ at position x and time t, C0 is the fixed concentration at x = 0, L is the length of the chamber, and D is the diffusion constant. Concentration vs. Distance and Time  0.8  C0  φ x,t  Normalized Concentration ( ( ) )  1.0  0.6  0.4  t=10.1s t=313.1s t=515.2s t=919.2s  0.2  0.0 0.0  0.2  0.4 0.6 0.8 Normalized Distance in Chamber ( ) x  1.0  L  Figure 3.21: Graphical depiction of the solution to the concentration profile in a one dimensional box of constant concentration reagent inflow as shown in Equation 3.4. Solved for a diffusion constant of 3 × 10−11 m2 s−1 with n = 100 terms. Imaging experiments were performed with inflow and outflow of ROX fluorescent dye (Invitrogen), monitoring the fluorescent intensity of cham60  3.3. Methods & Assay Development bers over time and visualizing concentration gradients across the length of the array to test the accuracy of the model (See Figure 3.22 for an example image of a device with a gradient due to insufficient reagent flow time). Results indicated the model does provide accurate estimates ±5 min for near uniform concentration of reagent in the chamber with inlet pressures of at least 9 psi and no downstream occlusion. Since removing the device from the microscope set-up for dye visualization is not trivial, a 2x safety margin is included in experiment design and the RT (and subsequent PCR) reagent flush is performed for 30 min.  Figure 3.22: Image of ROX fluorescent dye loaded device with insufficient flow time to achieve uniform reagent concentration across the array. Dye was loaded top to bottom at 7 psi for 10 min. The reverse transcription is performed for 55 °C for 1 h on a flat-bed thermal cycler. Comparison of various temperature RT reactions (See Section 3.4.1 showed that higher temperatures were more efficient, breaking the threshold in qPCR validation 5 cycles earlier than lower temperature (42 °C) RT (An enhancement factor of 32x). SuperScript™ III (Invitrogen) is one of the few reverse transcriptase enzymes capable of functioning at this tem61  3.3. Methods & Assay Development perature, with a measured activity half life of greater than 20 min at 55 °C (In comparison, SuperScript™ II, and M-MLV degrade to zero measurable activity after 5 min at 55 °C) [76]. In addition, SuperScript™ III has several engineered mutations to eliminate RNase H activity, increasing cDNA yield. For these reasons, it was chosen as the enzyme for reverse transcription. Polymerase Chain Reaction Following reverse transcription, the cDNA of interest is now affixed to the surface of the bead. The next step in a final design would be to remove the beads from the device for FACS of single beads into individual wells of a micro-well plate followed by multiplexed PCR for the heavy and light chains; this is easily accomplished by inverting the device and flowing out into a pipette tip. However, for purposes of device validation, the PCR step was performed directly on-chip to rapidly iterate through experiments, using probe or intercalating dye and a prototype BioMark™ (Fluidigm Corporation, San Francisco) system for qPCR imaging capabilities. The temperature of the device was returned to room temperature while PCR solution was prepared. For each primer set, the system was first tested with a TaqMan® probe to ensure specificity of PCR product on-device. A TaqMan® probe is an oligonucleotide selected to bind within the RNA region being amplified; on one end is affixed a fluorophore and the other end is linked to a quenching molecule. In this state fluorescence is suppressed. As PCR proceeds and DNA polymerase reaches the probe, the exonuclease activity cleaves the oligo, releasing the fluorophore from the quencher. The probe provides an added level of specificity as not only do the primers have to successfully bind to the template, but to see a signal resulting from PCR, the probe also has to bind to become cleaved. Following successful probe validation, EvaGreen® intercalating dye was used due to its lower cost. TaqMan® gene expression master mix (Invitrogen) was used to perform the on-device PCR. The PCR reaction mix with probe is shown in Table 3.2 (with half volume being used for in-tube reactions) and the reaction mix with intercalating dye is shown in Table 3.3.  62  3.3. Methods & Assay Development Reagent RNase-free Water TaqMan® Master Mix (2x) Primer Mix (20µM) Probe (20µM) Tween 1% BSA 1% Total Volume  Volume (µl) 13.38 25 1.0 0.625 5 5 50  Table 3.2: PCR reaction mixture with TaqMan® probe, used both for intube reactions (at half volume) and on-device (as shown). Reagent RNase-free Water Taqman® Master Mix (2x) Primer Mix (20µM) EvaGreen® (20x) Tween 1% BSA 1% Total Volume  Volume (µl) 13.38 25 1.0 2.500 5 5 50  Table 3.3: Polymerase chain reaction (PCR) reaction mixture with EvaGreen® , used both for in-tube reactions (at half volume) and on-device (as shown). PCR solution was loaded in the same way that RT solution is loaded; the device is flushed continuously for 30 min to ensure uniform distribution of reagent across the array. Following this, the secondary isolation valves are actuated to isolate the chambers. The device is then removed from the microscope while the valves are kept pressurized as it is transferred to the BioMark™ system; the system is equipped with a variable temperature vacuum chuck to maintain a constant position of the device relative to the lamp and lens apparatus while performing PCR. The device is then cycled through PCR with fluorescent imaging performed at the end of each annealing cycle.  63  3.4. Bead Testing Analysis The data resulting from on-chip RT-qPCR is a series of two images: one of the passive reference dye ROX contained in the Taqman® gene expression master mix, and one of the fluorescent dye or probe. Custom image analysis software was used to separate the image into individual chambers and assess the intensity change over time to build a qPCR curve for each chamber. The ROX dye is used as a baseline from which to normalize for intensity changes due to environmental and technical variability throughout the reaction. LabView Interfacing Custom LabView software was written to control solenoids for valve actuation and microscope control for imaging. The details are available upon request. The main interfaces used can be found in Appendix C.1.  3.4  Bead Testing  The bead-based approach was first validated in tube reactions followed by on-chip validation using quantitative PCR.  3.4.1  Validation in Tube Reactions with qPCR  Prior to on-device optimization, the bead-based approach was validated in tube reactions at equivalent concentrations to those expected in a working device. Endogenous control genes were chosen to test the protocol against simple transcripts and where probes could be designed to ensure the reaction product was specific. The genes GAPDH and b2m were chosen for initial testing on BaF3 cells (immortalized murine bone marrow B cells) for their high copy number. Primer/probe sets to generate products comparable in length to murine heavy chain transcripts were designed. For testing reverse transcription, the cells were lysed in tubes at 0.5 cell nl−1 and incubated with beads at 10 beads nl−1 at 55 °C with reverse  64  3.4. Bead Testing transcription mix as described in Section 3.3.26 . Several washing steps were performed followed by qPCR with either probe or intercalating dye. The efficiency of the dual functionalized beads was compared with poly-dT beads. The qPCR curves are shown in Figure 3.23 with the results summarized in Table 3.4.  Figure 3.23: Initial tube validation of dual-functionalized bead approach. Graph shows qPCR results for the various parameters described in Table 3.4. Results indicate that dual-functionalized beads work for successful reverse transcription and amplification of the target transcripts and are as efficient as commercially available poly-dT beads, though 32x less efficient than free primers in solution, possibly owing to steric effects or the spatial localization of available primer. The 55 °C reverse transcription incubation temperature was chosen after several optimization experiments with various temperatures and reverse transciptase enzymes were tried. Figure 3.24 shows the results of one experiment with SuperScript III. Results indicate a 55 °C RT reaction is approximately 32x more efficient than a 42 °C reaction or 50 °C reaction7 . It is hypothesized that this is due to the higher temperature helping the reverse transcriptase break through any RNA secondary structure. It is also noted that heat lysis reduces the 6 Beads were 50% covered with the appropriate primer and 50% covered with biotinylated BSA for testing purposes 7 Each cycle indicates a two-fold difference in concentration; therefore five cycles indicates a 25 = 32x efficiency difference  65  3.4. Bead Testing Description 10 beads nl−1 + EvaGreen read-out NTC for GAPDH 10 beads nl−1 + EvaGreen read-out NTC for b2m 10 beads nl−1 + EvaGreen read-out for GAPDH 10 beads nl−1 + EvaGreen read-out for b2m 10 beads nl−1 + EvaGreen read-out no 65 °C incubation and 50 °C RT for b2m 10 beads nl−1 + EvaGreen read-out no 65 °C incubation for GAPDH 10 beads nl−1 + EvaGreen read-out no 65 °C incubation for b2m 10 beads nl−1 + EvaGreen read-out poly-dT beads for GAPDH 10 beads nl−1 + EvaGreen readout no 65 °C incubation and 50 °C RT for GAPDH 10 beads nl−1 + Probe read-out NTC for b2m 10 beads nl−1 + Probe read-out for b2m 10 beads nl−1 + Probe read-out poly-dT beads for b2m 10 beads nl−1 + Probe read-out poly-dT beads for GAPDH Positive control EvaGreen read-out for GAPDH Positive control EvaGreen read-out for b2m Positive control Probe read-out for b2m  Ct NA NA 21.50 23.42 24.07 23.26 23.73 24.58 21.19 37.80 27.85 27.21 23.07 16.01 17.99 21.66  Table 3.4: Conditions and cycle threshold values for the various validation reactions carried out in tubes to support the validity of the dualfunctionalized bead approach. reaction efficiency by 16x in this experiment. This was also found to be the case on-chip as well; although not thoroughly investigated, it is thought that the high temperature denatures antibody on the bead surface, blocking otherwise active primer sites due to imaging experiments showing a loss of functional activity. The beads were also evaluated for their efficiency with varying proportional protein surface coverage. To enable adequate brightness for imaging and FACS, and to decrease the incubation time required for effective antibody capture, it is desired for the beads to have as many anti-IgG antibodies 66  3.4. Bead Testing  Figure 3.24: Comparison of different reverse transcription temperatures and heat lysis with bead’s surface functionalized with 50% biotinylated goat antimouse IgG and 50% b2m reverse primer. The cycle threshold (Ct ) values are shown in Table 3.5. Reaction Condition 55C RT with no heat lysis 42C RT 50C RT 55C RT with heat lysis No template control  Curve Colour Red Green Blue Orange Violet  Cycle Threshold Value 22.76 27.43 28.05 26.11 NA  Table 3.5: Ct values for Figure 3.24. on the surface as possible without significantly inhibiting the RT reaction or reducing the RT-PCR efficiency. Beads were functionalized with reverse strand primer for b2m and biotinylated goat anti-mouse IgG in 20% coverage increments and tested against single cell equivalent concentrations. The results are shown in Figure 3.25. The protein surface coverage can be increased to 60% before reaction efficiency starts noticeably decrease for a given incubation time. Note that this is the molar ratio loaded during bead functionalization. Due to the differing kinetics of binding between primer and protein, the realized surface ratio may be different than that in solution. These values thus essentially represent protein to mRNA ratios. Following successful validation of the RT-PCR reaction, it was necessary to validate the ability of the dual functionalized beads to retain the ability to bind target antibodies throughout the incubation and reverse transcription  67  3.4. Bead Testing  Cycle Threshold Value vs. Protein Coverage 10x Dilution 1x Dilution  Cycle Threshold Value  32  30  28  26  24 0  20  40 60 Protein Coverage (%)  80  100  Figure 3.25: Effect of proportional protein surface coverage on the RT-PCR efficiency of in-tube reactions at equivalent concentrations to device. The 1x dilution points are at a bead concentration of 10 beads nl−1 while the 10x dilution points are at a bead concentration of 1 beads nl. The 100% protein coverage bead samples come up late due to small amounts of genomic carryover between washing steps. steps.  68  3.4. Bead Testing  3.4.2  Antibody Binding & Imaging  To ensure antibody function remained intact throughout the binding and reaction steps and that it was sufficient for imaging and eventual FACS purposes, various binding validation experiments were performed. Beads of 20%, 40%, and 80% coverage of either goat anti-mouse immunoglobulin or lysozyme (the remaining covered with primer) were co-incubated with HyHel-5 murine hybridoma cells producing anti hen egg lysozyme (antiHEL) antibodies. Following incubation, the beads with goat anti-mouse IgG were co-incubated with fluorescently labelled HEL and the beads with bound lysozyme on the surface8 were co-incubated with labelled goat antimouse IgG. The results are shown in Table 3.6. It is observed that with all other variables being equal, the lysozyme covered beads, forming a direct assay for specificity, are several times more effective at generating signal than the IgG covered beads. This is likely because a labelled anti-mouse IgG would have multiple binding sites against the Fc region of the anti-HEL antibody bound to the surface lysozyme, whereas the indirect assay provides only two binding sites for lysozyme per antibody, one on each ‘arm’; for this reason, the lysozyme antigen covered beads are preferred for testing. However, it should be noted that for medium to low affinity antibodies, the antibody may be lost prior to sorting via FACS if antigen covered beads are used instead of high-affinity anti-IgG. There may be an additional avidity effect as well because the antibodies can interact with multiple surface-bound antigens. These considerations were not tested in this thesis. Following validation of binding and RT-PCR efficiency, the beads were also tested for binding capability following heat lysis. Heat lysis is the most elegant and simple way of lysing cells on device not needing the additional valving and flow structures required for chemical lysis. However, as shown in Table 3.7, heat lysis denatures the antibodies in solution to such a degree that binding is no longer effective. This is why chemical lysis was chosen for the protocol as described in Section 3.3.2 and the button valves designed. 8  Thus capturing the ‘arms’ of the HyHel-5 antibodies and exposing the Fc region for label binding  69  3.4. Bead Testing  Table 3.6: Antibody binding assay with varying protein surface coverage and either direct or indirect capture of specific antibodies. Bars represent measured fluorescent intensity.  Table 3.7: Assay for bead binding ability with 100% surface coverage of lysozyme against anti-HEL with and without a heat lysis step prior to reverse transcription. Subsequently, the detection limit of the beads was investigated to determine an estimate of the required incubation time prior to cell lysis. The experiment was set-up as described in Table 3.8 with results shown in Table 3.9. Detection was possible down to approximately 400 pM. A typical antibody secreting hybridoma or primary B cell can produce >300 antibodies per second; this gives an incubation time of approximately 27 min for a 2 nl chamber, matching the modelling results. Note that the negative controls of no surface coverage of lysozyme and a labelled anti-mouse antibody against lysozyme beads produced no non-specific absorption of labelled antibody. An example bead image is shown in Figure 3.26  70  3.4. Bead Testing  Table 3.8: Testing detection limit for functionalized beads against labelled antibodies. Results are shown in Table 3.9.  Tube Number A1 A2 A3 A5 A7 A9 C1 C2 C3  Intensity Above Background 70 84 38 53 63 0 0 77 0  Table 3.9: Results for the antibody binding reaction conditions shown in Figure 3.8  71  3.4. Bead Testing  Figure 3.26: Example 350 ms exposure of concentrated beads originally incubated at 1 bead per nL with labelled antibody at 4 nM concentration showing the ease of distinguishing signal from background.  72  3.5. Validation on Microfluidic Devices with qPCR  3.5  Validation on Microfluidic Devices with qPCR  Validation of the approach on microfluidic devices designed consisted of the approach outlined in Figure 3.17 and described in Section 3.3.2. For on-chip validation, a degenerate primer mix for murine IgG was used to synthesize cDNA for both heavy and light chains from HyHel-5 cells. Dual functionalized beads with the degenerate primer mix and goat anti-mouse IgG are loaded onto the microfluidic array after stochastic loading of cells9 . The cells are co-incubated with the beads following chemical lysis, and subsequently reverse transcription is performed to generate cDNA on the bead surface. Following reverse transcription, qPCR is performed to assay for successful cDNA synthesis. For controls, the individually addressible sub-arrays were used to perform parallel experiments. Sub-array one was loaded with both cells and beads, sub-array two with only beads, sub-array three with only cells, and sub-array four with nothing. As the chambers are isolated from each other, only those chambers containing both beads and cells should increase in fluorescence as PCR proceeds. The device is imaged during each annealing step of PCR. Image analysis is used to track the intensity of each chamber over time. This process yields a qPCR curve for each chamber as shown in Figure 3.27. In order to correlate the presence and number of cells with the measured qPCR curves, the contents of each chamber are manually counted and used to generate a heat-map of occupancy for the device as shown in Figure 3.28. A heat-map of qPCR values is also generated from the imaged intensity curves (Figure 3.29). These two heat-maps are superimposed to generate a data-set linking the occupancy with the cycle threshold value as shown in Figure 3.30. This data can then be used to generate box plots of cycle threshold value vs. occupancy number of each chamber in the array, as shown in Figure 3.31. As seen here, the cycle threshold number correlates with the occupancy of 9  Creates a poisson loading distribution  73  3.5. Validation on Microfluidic Devices with qPCR  Figure 3.27: Individual chamber qPCR curves from on-device PCR with imaging at each annealing step. The chosen threshold value is shown in red. the chambers as would be expected, with noise present due to inherent biological heterogeneity in antibody production and technical variability. Note that only 4 chambers out of the 704 chamber array with 0 cells came up at all. The results match the hypothesis for the experiment. Only those chambers which contained both cells and beads resulted in increasing fluorescent intensity from the qPCR validation step (except for 4 chambers which were nearby cell and bead containing chambers). The cycle threshold value also decreases with an increasing number of cells in a chamber as expected due to the increased number of transcript copies present. To further confirm the results observed via qPCR and ensure the product is specific, beads with mixed reverse primer for heavy and light chain murine antibodies were co-incubated with HyHel-5 hybridoma cells followed by chemical lysis and reverse transcription in tubes in parallel with a second on-device experiment. The beads on-device were extracted following reverse transcription and added to a tube containing PCR mix and primers 74  3.5. Validation on Microfluidic Devices with qPCR  Figure 3.28: Example of manually curated heat-map for cell-count occupancy of each chamber in the array. The lanes were loaded vertically. In this experiment, subarray 1 (far left) was a positive control test with both beads and cells, subarray 2 and 4 were empty, and subarray 3 was a negative control test with cells but no beads (and therefore no RT primer). for nested PCR in parallel with the in-tube reactions. The products were added to an agarose gel for gel electrophoresis; the resulting image is shown in Figure 3.32. Both heavy and light chain bands are observed indicating a successful synthesis of cDNA on the surface of the bead. Thus the prototype device designed is shown to be a successful approach to single-cell bead-based capture of both sequence and antibody specificity information. Future work should investigate scaling up the design and incorporating downstream sorting via FACS and subsequent cloning and recombinant expression followed by detailed kinetics measurements.  75  3.5. Validation on Microfluidic Devices with qPCR  Figure 3.29: qPCR cycle threshold value heat map for the device with occupancy shown in Figure 3.28. As expected, only subarray 1 is coming up during PCR.  76  3.5. Validation on Microfluidic Devices with qPCR  Figure 3.30: Superimposed heat-map of cell occupancy count (white circles, from Figure 3.28) with the qPCR cycle threshold information of Figure 3.29.  77  3.5. Validation on Microfluidic Devices with qPCR  Figure 3.31: Box plot of cycle threshold value vs. occupancy number of each chamber in the array.  78  3.5. Validation on Microfluidic Devices with qPCR  Figure 3.32: Agarose gel electrophoresis of RT-PCR products. Lane 1 through lane 5 contain beads with protein surface coverage of 0%, 20%, 40%, 60%, and 80% respectively, with RT-PCR performed in tubes. Lane 6 contains beads with 60% primer coverage with reverse transcription performed on-chip. Lane 7 contains a 1 cell nl−1 equivalent concentration of RNA against free primers, and lane 8 contains a 0.1 cell nl−1 equivalent concentration of RNA against free primers. Lane 9 contains beads from the device that were not exposed to cells, and lane 10 contains beads from tubes that were not exposed to cells.  79  Chapter 4  Conclusions & Recommendations The goals of this thesis were: to design, develop and validate a microfluidic device and assay for pre-screening for targeted antibodies from single antibody secreting cells; to ensure the device was easily scalable, robust, and amenable to automation; and to design and code a numerical modelling pipeline to enable the characterization and optimization of microfluidic designs for solid substrate capture. The first two goals were achieved by designing and validating an assay protocol and microfluidic device allowing for parallel stochastic loading of single antibody secreting cells into a planar array of micro-wells. A novel design incorporating a large planar array and button valve architecture was created and validated. Co-incubating with dual-functionalized beads with reverse transcription of transcript mRNA allowed for simultaneous capture of affinity and heavy and light chain sequence information, verified with qPCR and gel electrophoresis, allowing for downstream sorting and sequencing. An off-chip and on-chip protocol was created and optimized for ensuring a robust and repeatable pipeline. The device architecture utilized computer controlled integrated pneumatic valves, and the repeating unit-cell structure of the array lends itself to automated chamber tracking via image analysis software. The goal for a numerical modelling pipeline was achieved by the coding and development of a MATLABCOMSOL integrated pipeline to allow for user-customized investigation into the effect of chamber geometry, incubation parameters, and bead distributions in diffusion-based solid substrate capture of secreted molecules. The system is robust to user input and allows for extensive characterization prior  80  Chapter 4. Conclusions & Recommendations to any device fabrication. The modelling was able to accurately predict incubation periods and bead flux distributions in microfluidic chambers. The results indicate that future work should investigate more accurately controlling the bead quantity and distribution within the chambers to maximize homogeneity of surface protein capture. While enabling high-throughput screening and sequence information coupling with scaling, the approach developed in this thesis only serves to function as a pre-screen and cannot yet replace existing downstream approaches for validating kinetics and cross-reactivity. As a large sample size of antibody secreting cells will have a wide range of affinities to a target antigen [77], the resulting bead fluorescent intensity will have a wide continuous range. In FACS, a cut-off threshold will be set to retain beads above an arbitrary intensity [78, 79]. However, absent any additional information, this is a relative measurement between beads and is only correlated with affinity via some scaling factor based on incubation time, bead distribution, chamber geometry, secretion rate, and experimental noise. This device and approach would need to be coupled with downstream screening approaches such as microfluidic kinetics screening, immunoprecipitation, and/or ELISA to test for high affinity and cross-reactivity [14, 80–82]. Future work on this design could investigate the possibility of using multiple antigens labelled with spectrally distinct fluorophores to allow for multiplexed FACS elimination of antibodies displaying cross-reactivity. Recent work has shown the ability to screen single antibody secreting cells for high affinity antibodies while simultaneously obtaining kinetic information, all in an automated microfluidic device [14]. While a breakthrough technology, this approach requires time-consuming serial analysis of each individual cell, prior to any knowledge of affinity. An approach integrating bead-based assays with planar parallel arrays for rapid stochastic loading is also under development [14, 54]; while effective, this approach is still limited in throughput at this time to thousands of micro-wells. The approaches developed to date provide excellent methods for screening for high-affinity antibodies from an enriched sample population but lack scalability to large sample sizes and parallel operation to pre-screen relevant antibodies out of 81  Chapter 4. Conclusions & Recommendations a pool of unknowns [48–52]. The micro-well approaches developed thus far require manual pipetting to either load, or to extract the antibody secreting cells [48, 49] limiting the achievable throughput. The emulsion linkage-PCR approach of DeKosky et. al. [52] has the potential for extremely highthroughput sequencing of the antibody repertoire of a given sample but still requires that all genes be cloned and expressed prior to affinity screening, drastically increasing the workload for large sample sizes, similar to the hybridoma method [14–16, 18, 52]. We note that this approach presented herein provides a means to solve this issue; it may be combined with next generation sequencing of antibody repertoires to allow for the assignment of binding specificity to different antibody sequences and sequence clusters.  82  Works Cited [1] J. Playfair and G. Bancroft, infection & immunity, 4th ed.  Great  Clarendon Street, Oxford, OX2 6DP, United Kingdom: Oxford University Press, 2013. [2] K. M. Murphy, P. Travers, and M. Walport, Janeway’s Immunobiology.  Garland Science / Taylor & Francis LLC, 2011.  [3] P. J. Carter, “Potent antibody therapeutics by design,” Nature Reviews Immunology, vol. 6, 2006. [Online]. Available:  http:  //dx.doi.org/10.1038/nri1837 [4] T. Dharakul, S. Songsivilai, N. Anuntagool, W. Chaowagul, S. Wongbunnate, P. Intachote, and S. Sirisinha, “Diagnostic value of an antibody enzyme-linked immunosorbent assay using affinity-purified antigen in an area endemic for melioidosis,” American Journal of Tropical Medicine and Hygiene, vol. 56, no. 4, pp. 418–423, 1997. [5] C. L. Reading and Y. Takaue, “Monoclonal antibody applications in bone marrow transplantation,” Biochimica et Biophysica Acta, vol. 865, no. 2, pp. 141–170, October 1986. [6] S. F. Kingsmore, “Multiplexed protein measurement: technologies and applications of protein and antibody arrays,” Nature Reviews Drug Discovery, vol. 5, pp. 310–321, April 2006. [Online]. Available: http://dx.doi.org/10.1038/nrd2006 [7] D. J. Cahill, “Protein and antibody arrays and their medical applications,” Journal of Immunological Methods, vol. 250, no. 1-2, pp. 81–91, 2001. [Online]. Available:  http://dx.doi.org/10.1016/  S0022-1759(01)00325-8 [8] B. L. Ackermann and M. J. Berna, “Coupling immunoaffinity techniques with MS for quantitative analysis of low-abundance protein 83  Works Cited biomarkers,” Expert Reviews of Proteomics, vol. 4, pp. 501–538, 2007. [9] A. Bradbury, N. Velappan, V. Verzillo, M. Ovecka, L. Chasteen, D. Sblattero, R. Marzari, J. Lou, R. Siegel, and P. Pavlik, “Antibodies in proteomics II: screening, high-throughput characterization and downstream applications,” Trends Biotechnology, vol. 21, no. 7, pp. 312–317, July 2003. [10] J. C. Love, “Making antibodies from scratch,” Nature Biotechnology, vol. 28, pp. 1176–1178, 2010. [Online]. Available: http://dx.doi.org/ 10.1038/nbt1110-1176 [11] M. Little, “Antibody Libraries: How to Overcome the Discovery LogJam,” Biotechnology, 2006. [12] M. Clark, “Antibody humanization: a case of the ‘Emporer’s new clothes’,” Immunology Today, vol. 21, no. 8, pp. 397–402, August 2000. [Online]. Available: http://dx.doi.org/10.1016/S0167-5699(00) 01680-7 [13] T. Tiller, C. E. Busse, and H. Wardemann, “Cloning and expression of murine Ig genes from single B cells,” Journal of Immunological Methods, vol. 350, pp. 183–193, August 2009. [Online]. Available: http://dx.doi.org/10.1016/j.jim.2009.08.009 [14] A. Singhal, “Microfluidic Technologies for Rapid, High-Throughput Screening and Selection of Monoclonal Antibodies from Single Cells,” Ph.D. dissertation, Univ. of British Columbia, Vancouver, British Columbia, November 2012. [15] G. K¨ ohler and C. Milstein, “Continuous cultures of fused cells secreting antibody of predefined specificity,” Nature, vol. 256, pp. 495–497, August 1975. [Online]. Available: http://dx.doi.org/10.1038/256495a0 [16] C. Milstein, “The hybridoma revolution:  an offshoot of basic  research,” BioEssays, vol. 21, no. 11, pp. 966–973, November 1999. [Online]. Available:  http://dx.doi.org/10.1002/(SICI)  1521-1878(199911)21:11 966::AID-BIES9 3.0.CO;2-Z [17] C. L. Cepko, P. S. Changelian, and P. A. Sharp, “Immunoprecipitation with two-dimensional pools as a hybridoma screening technique: Production and characterization of monoclonal antibodies against 84  Works Cited adenovirus 2 proteins,” Virology, vol. 110, no. 2, pp. 385–401, April 1981. [Online]. Available: http://dx.doi.org/10.1016/0042-6822(81) 90069-6 [18] C. Cervino, E. Weber, D. Knopp, and R. Niessner, “Comparison of hybridoma screening methods for the efficient detection of high-affinity hapten-specific monoclonal antibodies,” Journal of Immunological Methods, vol. 329, no. 1-2, pp. 184–193, January 2008. [Online]. Available: http://dx.doi.org/10.1016/j.jim.2007.10.010 [19] I. Sogut, I. Hatipoglu, G. Kanbak, and A. Basalp, “Monoclonal antibodies specific for hepatitis B e antigen and hepatitis B core antigen,” Hybridoma, vol. 30, no. 5, pp. 475–479, October 2011. [Online]. Available: http://dx.doi.org/10.1089/hyb.2011.0052 [20] G. M. Whitesides, “The origins and future of microfluidics,” Nature, vol. 442, pp. 368–373, 2006. [Online]. Available:  http:  //dx.doi.org/10.1038/nature05058 [21] A. Manz, D. J. Harrison, E. M. Verpoorte, J. C. Fettinger, A. Paulus, H. L¨ udi, and H. M. Widmer, “Planar chips technology for miniaturization and integration of separation techniques into monitoring systems: Capillary electrophoresis on a chip,” Journal of Chromatography A, vol. 593, no. 1-2, pp. 253–258, February 1992. [Online]. Available: http://dx.doi.org/10.1016/0021-9673(92)80293-4 [22] E. Bassous, H. H. Taub, and L. Kuhn, “Ink jet printing nozzle arrays etched in silicon,” Applied Physics Letters, vol. 31, 1977. [Online]. Available: http://dx.doi.org/10.1063/1.89587 [23] M. J. Zdeblick, P. W. Barth, and J. B. Angell, “Microminiature fluidic amplifier,” Social Studies of Science, 1986. [24] K. A. Hagan,  W. L. Meier,  J. P. Ferrance,  and J. P.  Landers, “Chitosan-Coated Silica as a Solid Phase for RNA Purification in a Microfluidic Device,”  Analytical Chemistry,  vol. 81, no. 13, pp. 5249–5256, June 2009. [Online]. Available: http://dx.doi.org/10.1021/ac900820z [25] W.-W. Jeong and C. Kim, “One-step method for monodisperse microbiogels by glass capillary microfluidics,” Colloids and Surfaces 85  Works Cited A: Physicochemical and Engineering Aspects, vol. 384, no. 1-3, pp. 268–273, July 2011. [Online]. Available: http://dx.doi.org/10.1016/j. colsurfa.2011.04.006 [26] C.-H. Lin, G.-B. Lee, B.-W. Chang, and G.-L. Chang, “A new fabrication process for ultra-thick microfluidic microstructures utilizing SU-8 photoresist,” Journal of Micromechanics and Microengineering, vol. 12, no. 5, p. 590, 2002. [Online]. Available: http://dx.doi.org/10.1088/0960-1317/12/5/312 [27] Y. Xia and G. M. Whitesides, “Soft Lithography,” Annual Review of Materials Science, vol. 28, pp. 153–184, 1998. [28] M. A. Unger, H.-P. Chou, T. Thorsen, A. Scherer, and S. R. Quake, “Monolithic Microfabrabricated Valves and Pumps by Multilayer Soft Lithography,” Science, vol. 288, no. 5463, pp. 113–116, 2000. [Online]. Available: http://dx.doi.org/10.1126/science.288.5463.113 [29] J. C. McDonald and G. M. Whitesides, “Poly(dimethysiloxane) as a Material for Fabricating Microfluidic Devices,” Accounts of Chemical Research, vol. 35, no. 7, July 2002. [30] T. M. Squires and S. R. Quake, “Microfluidics: fluid physics at the nanolitre scale,” Reviews of Modern Physics, vol. 77, pp. 977–1026, October 2005. [31] J. C. McDonald, D. C. Duffy, J. R. Anderson, D. T. Chiu, H. Wu, O. J. A. Schueller, and G. M. Whitesides, “Fabrication of microfluidic systems in poly(dimethylsiloxane),” Electrophoresis, vol. 21, pp. 27– 40, 2000. [32] Joshua S. Marcus, W. French Anderson, and Stephen R. Quake, “Microfluidic Single-Cell mRNA Isolution and Analysis,” Analytical Chemistry, vol. 78, no. 9, pp. 3084–3089, 2006. [Online]. Available: http://pubs.acs.org/doi/full/10.1021/ac0519460 [33] A. R. Wheeler,  W. R. Throndset,  R. J. Whelan,  A. M.  Leach, R. N. Zare, Y. H. Liao, K. Farrell, I. D. Manger, and A. Daridon, “Microfluidic Device for Single-Cell Analysis,” Analytical Chemistry, vol. 75, no. 14, pp. 3581–3586, 2003. [Online]. Available: http://dx.doi.org/10.1021/ac0340758 86  Works Cited [34] A. K. White,  M. VanInsberghe,  O. I. Petriv,  M. Hamidi,  D. Sikorski, M. A. Marra, J. Piret, S. Aparicio, and C. L. Hansen, “High-throughput microfluidic single-cell RT-qPCR,” PNAS, vol. 108, no. 34, pp. 13 999–14 004, August 2010. [Online]. Available: http://dx.doi.org/10.1073/pnas.1019446108 [35] K. Leung, H. Zahn, T. Leaver, K. M. Konwar, N. W. Hanson, A. P. Pag´e, C.-C. Lo, P. S. Chain, S. J. Hallam, and C. L. Hansen, “A programmable droplet-based microfluidic device applied to multi parameter analysis of single microbes and microbial communities,” PNAS, vol. 109, no. 20, pp. 7665–7670, February 2012. [Online]. Available: http://dx.doi.org/10.1073/pnas.1106752109 [36] A. Singhal, C. A. Haynes, and C. L. Hansen, “Microfluidic Measurements of Antibody-Antigen Binding Kinetics from LowAbundance Samples and Single Cells,” Analytical Chemistry, vol. 80, no. 20,  pp. 8671–8679,  September 2010. [Online]. Available:  http://dx.doi.org/10.1021/ac101956e [37] J. J. Agresti, E. Antipov, A. R. Abate, K. Ahn, A. C. Rowat, J.-C. Baret, M. Marquez, A. M. Mlibanov, A. D. Griffiths, and D. A. Weitz, “Ultrahigh-throughput screening in drop-based microfluidics for directed evolution,” PNAS, vol. 107, no. 9, March 2010. [Online]. Available: http://dx.doi.org/10.1073/pnas.0910781107 [38] S. Bhattacharya, A. Datta, J. M. Berg, and S. Gangopadhyay, “Studies on Surface Wettability of Poly(Dimethyl) Siloxane (PDMS) and Glass Under Oxygen-Plasma Treatment and Correlation With Bond Strength,” Journal of Microelectromechanical Systems, vol. 14, no. 3, pp. 590–597, June 2005. [39] M. A. Eddings, M. A. Johnson, and B. K. Gale, “Determing the optimal PDMS-PDMS bonding technique for microfluidic devices,” Journal of Micromechanics and Microengineering, vol. 18, no. 6, 2008. [Online]. Available: http://dx.doi.org/10.1088/0960-1317/18/ 6/067001 [40] B.-H. Jo, L. M. V. Lerberghe, K. M. Motsegood, and D. J. Beebe, “Three-dimensional micro-channel fabrication in polydimethylsiloxane 87  Works Cited (pdms) elastomer,” Journal of Microelectromechanical Systems, vol. 9, no. 1, March 2000. [41] D. Whitcombe, J. Theaker, S. Guy, T. Brown, and S. Little, “Detection of PCR products using self-probing amplicons and fluorescence,” Nature Biotechnology, vol. 17, no. 8, pp. 804–807, August 1999. [42] S. Tyagi and F. Kramer, “Molecular beacons: probes that fluoresce upon hybridization,” Nature Biotechnology, vol. 14, no. 3, pp. 303–308, March 1996. [43] T. Morrison, J. Weis, and C. Wittwer, “Quantification of low-copy transcripts by continuous SYBR Green I monitoring during amplification,” Biotechniques, vol. 24, no. 6, pp. 954–958, June 1998. [44] R. Saiki, D. Gelfand, S. Stoffel, S. Scharf, R. Higuchi, G. Horn, K. Mullis, and H. Erlich, “Primer-directed enzymatic amplification of DNA with a thermostable DNA polymerase,” Science, vol. 239, no. 4839, pp. 487–491, January 1988. [45] J. Wrammert, K. Smith, J. Miller, W. A. Langley, K. Kokko, C. Larsen, N.-Y. Zheng, I. Mays, L. Garman, C. Helms, J. James, G. M. Air,  J. D. Capra,  R. Ahmed,  and P. C. Wilson,  “Rapid cloning of high-affinity human monoclonal antibodies against influenza virus,” Nature, vol. 453, May 2008. [Online]. Available: http://dx.doi.org/10.1038/nature06890 [46] D. Corti, J. Voss, S. J. Gamblin, G. Codoni, A. Macagno, D. Jarrossay, S. G. Vachieri, D. Pinna, A. Minola, F. Vanzetta, C. Silacci, B. M. Fernandez-Rodriguez, G. Agatic, S. Bianchi, I. Giacchetto-Sasselli, L. Calder, F. Sallusto, P. Collins, L. F. Haire, N. Temperton, J. P. Skehel, and A. Lanzavecchia, “A neutralizing antibody selected from plasma cells that binds to group 1 and group 2 influenza A hemagglutinins,” Science, vol. 333, no. 6044, pp. 850–856, 2011. [Online]. Available: http://dx.doi.org/10.1126/science.1205669 [47] J. S. Babcook, K. B. Leslie, O. A. Olsen, R. A. Salmon, and J. W. Schrader, “A novel strategy for generating monoclonal antibodies from single, isolated lymphocytes producing antibodies of defined specificities,” PNAS, vol. 23, no. 93, pp. 7843–7848, July 1996. 88  Works Cited [48] J. C. Love, J. L. Ronan, G. M. Grotenbreg, A. G. van der Even, and H. L. Ploegh, “A microengraving method for rapid selection of single cells producing antigen-specific antibodies,” Nature Biotechnology, vol. 24, no. 6, pp. 703–707, June 2006. [49] A. Jin, T. Ozawa, K. Tajiri, T. Obata, K. Kinoshita, S. Kadawaki, K. Takahashi, T. Sugiyama, H. Kishi, and A. Muraguchi, “A rapid and efficient single-cell manipulation method for screening antigenspecific antibody-secreting cells from human peripheral blood,” Nature Medicine, vol. 15, no. 9, pp. 1088–1092, 2009. [Online]. Available: http://dx.doi.org/10.1038/nm.1966 [50] C. M. Story, E. Papa, C.-C. A. Hu, J. L. Ronan, K. Herlihy, H. L. Ploegh, and J. C. Love, “Profiling antibody responses by multi parametric analysis of primary b cells,” PNAS, vol. 105, no. 46, November 2009. [Online]. Available: http://dx.doi.org/10.1073/pnas. 0805470105 [51] S. K¨ oster, F. E. Angil`e, H. Duan, J. J. Agresti, A. Wintner, C. Schmitz, A. C. Rowat, C. A. Merten, D. Pisignano, A. D. Griffiths, and D. A. Weitz, “Drop-based microfluidic devices for encapsulation of single cells,” Lab on a Chip, vol. 8, pp. 1110–1115, May 2008. [Online]. Available: http://dx.doi.org/10.1039/B802941E [52] B. J. DeKosky, G. C. Ippolito, R. P. Deschner, J. J. Lavinder, Y. Wine, B. M. Rawlings, N. Varadarajan, C. Giesecke, T. D¨orner, S. F. Andrews, P. C. Wilson, S. P. Hunicke-Smith, C. G. Willson, A. D. Ellington, and G. Georgiou, “High-throughput sequencing of the paired human immunoglobulin heavy and light chain repertoire,” Nature Biotechnology, vol. 31, pp. 166–169, January 2013. [Online]. Available: http://dx.doi.org/10.1038/nbt.2492 [53] M. G. Carter, A. A. Sharov, V. VanBuren, D. B. Dudekula, C. E. Carmack, C. Nelson, and M. S. Ko, “Transcript copy number estimation using a mouse whole-genome oligonucleotide microarray,” Genome Biology, June 2005. [Online]. Available: http://dx.doi.org/10.1186/gb-2005-6-7-r61 [54] V. Lecault, M. VanInsberghe, S. Sekulovic, D. J. H. F. Knapp, 89  Works Cited S. Wohrer, W. Bowden, F. Viel, T. McLaughlin, A. Jarandehei, M. Miller, D. Falconnet, A. K. White, D. G. Kent, M. R. Copley, F. Taghipour, C. J. Eaves, K. Humphries, J. M. Piret, and C. L. Hansen, “High-throughput analysis of single hematopoietic stem cell proliferation in microfluidic cell culture arrays,” Nature Methods, vol. 8, pp. 581–586, May 2011. [Online]. Available: http://dx.doi.org/10.1038/nmeth.1614 [55] (2013) Multiphysics Modeling and Simulation Software. [Online]. Available: http://www.comsol.com/ [56] COMSOL Multiphysics Modeling Guide, PDF Manual, COMSOL, September 2005, v3.2. [57] COMSOL Multiphysics Modeling Guide, PDF Manual, COMSOL, October 2007, v3.4. [58] COMSOL Multiphysics User’s Guide, PDF Manual, COMSOL, October 2010, v4.1. [59] COMSOL Multiphysics Reference Guide, PDF Manual, COMSOL, May 2010, v4.3. [60] COMSOL Multiphysics User’s Guide, PDF Manual, COMSOL, May 2012, v4.3. [61] O. C. Ziekiewicz and R. L. Taylor, The Finite Element Method Set. Jordan Hill, GBR: Butterworth-Heinemann, 2005. [62] J. N. Reddy, Finite Element Method, 2nd ed.  New York: McGraw-  Hill, 1993. [63] O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method, 5th ed. Linacre House, Jordan Hill, Oxford: Butterworth-Heinemann, 2000. [64] W. B. J. Zimmerman, Process Modelling and Simulation with Finite Element Methods.  Singapore: World Scientific, 2004.  [65] G. R. Liu and S. S. Quek, Finite Element Method: A Practical Course. Jordan Hill, GBR: Butterworth-Heinemann, May 2003. [66] F.-J. Sayas, “A gentle introduction to the Finite Element Method,” 2008.  [Online].  Available:  ftp://tug.ctan.org/pub/tex-archive/  macros/latex/contrib/IEEEtran/bibtex/IEEEtran bst HOWTO.pdf 90  Works Cited [67] (2013) MATLAB — The Language of Technical Computing. [Online]. Available: http://www.mathworks.com/products/matlab/ [68] R. M. Robertson, S. Laib, and D. E. Smith, “Diffusion of isolated DNA molecules:  Dependence on length and topology,” PNAS,  vol. 103, no. 19, pp. 7310–7314, May 2006. [Online]. Available: http://dx.doi.org/10.1073/pnas.0601903103 [69] A. M. Yoffe, P. Prinson, A. Gopal, C. M. Knobler, W. M. Gelbart, and A. Ben-Shaul, “Predicting the sizes of large RNA molecules,” PNAS, vol. 105, no. 42, pp. 16 153–16 158, August 2008. [Online]. Available: http://dx.doi.org/10.1073/pnas.0808089105 [70] K. A. Heyries, C. Tropini, M. VanInsberghe, C. Doolin, O. I. Petriv, A. Singhal, K. Leung, C. B. Hughesman, and C. L. Hansen, “Megapixel digital PCR,” Nature Methods, vol. 8, pp. 649–651, 2011. [Online]. Available: http://dx.doi.org/10.1038/nmeth.1640 [71] H. Bruus, Theoretical Microfluidics. Great Clarendon Street, Oxford, OX2 6DP: Oxford University Press, 2008. [72] W. M. Haynes, D. R. Lide, and T. J. Bruno, CRC Handbook of Chemistry and Pphysics 2012-2013, 93rd ed. 6000 Broken Sound Parkway, NW, Suite 300, Boca Raton, FL 33487-2742: CRC Press, 2012. [73] G. Woodbury, Introduction to Statistics.  Cengage Learning, 2001.  [74] T. Yamada, H. Soma, and S. Morishita, “Primerstation: a highly specific multiplex genomic PCR primer design server for the human genome,” Nucleic Acids Research, vol. 34, pp. 665–669, July 2006. [Online]. Available: http://dx.doi.org/10.1093/nar/gkl297 [75] R. Owczarzy, B. G. Moreira, Y. You, M. A. Behlke, and J. A. Walder, “Predicting Stability of DNA Duplexes in Solutions Containing Magnesium and Monovalent Cations,” Biochemistry, vol. 47, no. 19, pp. 5336–5353, April 2008. [Online]. Available: http://dx.doi.org/10.1021/bi702363u [76] (2013) Available:  Superscript  III  Reverse  Transcriptase.  [Online].  http://www.invitrogen.com/site/us/en/home/  Products-and-Services/Applications/PCR/reverse-transcription/ reverse-transcriptase-enzymes/superscript-iii-reverse-transcriptase. 91  Works Cited html [77] F. Karush, “The Affinity of Antibody: Range, Variability, and the Role of Multivariance,” Immunoglobulins Comprehensive Immunology, vol. 5, pp. 85–116, 1978. [78] W. Bonner, H. Hulett, R. G. Sweet, and L. Herzenberg, “Fluorescence Activated Cell Sorting,” Review of Scientific Instruments, vol. 43, no. 3, 1972. [Online]. Available: http://dx.doi.org/10.1063/1.1685647 [79] J. Sandberg, B. Werne, M. Dessign, and J. Lundeberg, “Rapid flow-sorting to simultaneously resolve multiplex massively parallel sequencing products,” Scientific Reports, no. 108, October 2011. [Online]. Available: http://dx.doi.org/10.1038/srep00108 [80] M.  Selbach  and  M.  Mann,  “Protein  interaction  screening  by quantitative immunoprecipitation combined with knockdown (QUICK),” Nature Methods, pp. 981–983, 2006. [Online]. Available: http://dx.doi.org/10.1038/nmeth972 [81] R. M. T. de Wildt, C. R. Mundy, B. D. Gorick, and I. M. Tomlinson, “Antibody arrays for high-throughput screening of antibody-antigen interactions,” Nature Biotechnology, vol. 18, pp. 989–994, 2000. [Online]. Available: http://dx.doi.org/10.1038/79494 [82] J. C. Miller, H. Zhou, J. Kwekel, R. Cavallo, J. Burke, E. B. Butler, B. S. The, and B. B. Haab, “Antibody microarray profiling of human prostate cancer sera: Antibody screening and identification of potential biomarkers,” Proteomics, vol. 3, no. 1, 2003. [Online]. Available: http://dx.doi.org/10.1002/pmic.200390009 [83] A. S. Bedekar, Y. Wang, S. S. Siddhaye, S. Krishnamoorthy, and S. F. Malin, “Design Software for Application-Specific Microfluidic Devices,” Clinical Chemistry, vol. 53, no. 11, pp. 2023–2026, November 2007. [Online]. Available:  http://dx.doi.org/10.1373/  clinchem.2007.090498 [84] T. Thorsen, S. J. Maerkl, and S. R. Quake, “Microfluidic Large-Scale Integration,” Science, vol. 298, no. 5593, pp. 580–584, October 2002. [Online]. Available: http://dx.doi.org/10.1126/science.1076996 [85] Jessica Melin and Stephen R. Quake, “Microfluidic Large-Scale 92  Works Cited Integration:  The Evolution of Design Rules for Biological  automation,”  Annual Review of Biophysics and Biomolecular  Structure, vol. 36, pp. 213–231, 2007. [Online]. Available:  http:  //dx.doi.org/10.1146/annurev.biophys.36.040306.132646 [86] E. Leclerc, Y. Sakai, and T. Fujii, “A multi-layer PDMS microfluidic device for tissue engineering applications,” Transducers, pp. 415–418, 2003. [87] G. Vozzi, C. Flaim, A. Ahluwalia, and S. Bhatia, “Fabrication of PLGA scaffolds using soft lithography and microsyringe deposition,” Biomaterials, vol. 24, pp. 2533–2540, 2003. [88] M. Zhang, J. Wu, L. Wang, K. Xiao, and W. Wen, “A simple method for fabricating multi-layer PDMS structure for 3D microfluidic chips,” Lab on a Chip, vol. 10, pp. 1199–1203, February 2010. [89] J. Huft, D. J. D. Costa, D. Walker, and C. L. Hansen, “Tthreedimensional large-scale microfluidic integration by laser ablation of interlayer connections,” Lab on a Chip, vol. 10, pp. 2358–2365, June 2010. [90] M. W. Toepke, V. V. Abhyankar, and D. J. Beebe, “Microfluidic logic gates and timers,” Lab on a Chip, vol. 7, pp. 1449–1453, 2007. [Online]. Available: http://dx.doi.org/10.1039/B708764K [91] J. A. Weaver, J. Melin, D. Stark, S. R. Quake, and M. A. Horowitz, “Static control logic for microfluidic devices using pressure-gain valves,” Nature Physics, vol. 6, pp. 218–223, January 2010. [Online]. Available: http://dx.doi.org/10.1038/nphys1513 [92] J.-C. Baret, O. J. Miller, V. Taly, M. Ryckelynck, A. El-Harrak, L. Frenz, C. Rick, M. L. Samuels, J. B. Hutchison, J. J. Agresti, D. R. Link, D. A. Weitz, and A. D. Griffiths, “Fluorescent-activated droplet sorting (FADS): efficient microfluidic cell sorting based on enzymatic activity,” Lab on a Chip, April 2009. [93] C. Ricciardi, G. Canavese, R. Castagna, I. Ferrante, A. Ricci, S. L. Marasso, L. Napione, and F. Bussolino, “Integration of microfluidic and cantilever technology for biosensing application in liquid environment,” Biosensors and Bioelectronics, vol. 26, pp. 93  Works Cited 1565–1570, May 2010. [Online]. Available:  http://dx.doi.org/10.  1016/j.bios.2010.07.114 [94] A. Nadim, “Modeling of mass transfer limitation in biomolecular assays.” Annals of the New York Academy of Sciences, vol. 1161, pp. 34–43, 2009. [Online]. Available: http://dx.doi.org/10.1111/j. 1749-6632.2008.04071.x [95] G. James, D. Burley, D. Clements, P. Dyke, J. Searl, N. Steele, and J. Wright, Advanced Modern Engineering Mathematics, 3rd ed. Essex, England: Pearson, 2004. [96] M. Shapiro and H. Brenner, “Dispersion of a chemically reactive solute in a spatially periodic model of a porous medium,” Chemical Engineering Science, vol. 43, no. 3, pp. 551–571, 1988. [97] D. A. Edwards, M. Shapiro, and H. Brenner, “Dispersion and reaction in two-dimensional model porous media,” Physics of Fluids, vol. 5, no. 4, p. 837, 1993. [98] B. A. Finlayson and L. E. Scriven, “The Method of Weighted Residuals — A Review,” Applied Mechanics Reviews, vol. 19, no. 9, pp. 735–748, September 1966. [99] C. Greenough. (2001) The Finite Element Library - Theory and Programming Techniques - Introduction to Time-Dependent Problems. [Online]. Available: http://www.softeng.rl.ac.uk/st/projects/felib3/ Docs/html/Intro/intro-node45.html [100] J. D. Hoffman and S. Frankel, Numerical Methods for Engineers and Scientists, 2nd ed. 270 Madison Avenue, New York, NY 10016: CRC Press, May 2001. [101] P. Niyogi, Introduction to Computational Fluid Dynamics.  India:  Pearson Education, 2006. [102] D. P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods.  Academic Press, 1996.  [103] ——, Nonlinear Programming, 2nd ed.  Athena Scientific, 1999.  [104] (2012) Lagrange multipliers. Encyclopedia of Mathematics. [Online]. Available:  http://www.encyclopediaofmath.org/index.php?  title=Lagrange multipliers&oldid=24258 94  [105] C. Greenough. (2001) The Finite Element Library - Theory and Programming Techniques - Lagrange Multipliers. [Online].  Available:  http://www.softeng.rl.ac.uk/st/projects/felib3/  Docs/html/Intro/intro-node68.html  95  Appendix A  Numerical Modelling Pipeline One of the outstanding problems with microfluidic design is the distinct lack of software tools devoted to helping the interested engineer, chemist, or biologist design a device and predict its performance for a specified task [83]. Microfluidics layout and design has reached the point of complexity of early microelectronic circuits [84, 85] with multiple stacked interconnecting planar layers [86–88], complex interlayer connections by laser ablation [89], valving logic [84, 90, 91] and multiplexed arrays [84, 85], and even integrated electronic components [92, 93]. Despite this increasing level of sophistication, the layout software, physics, and performance modelling tools available to circuit designers have no analogue for microfluidics. While design rules and common practises exist [84, 85], much of the design work involves extensive experience and intuition, and often many iterations of prototypes are necessary to elucidate device parameters and achieve a functional design; this design must often then be further optimized to achieve the necessary performance. During the course of design investigation, several problems arose whose solutions were not easily answered experimentally. The ultimate driving motivation for this modelling is to accurately calculate and keep track of flux on a given boundary in a microfluidic chamber over time. The opportunity was used to create a mathematical approach, software system, and associated code that, while limited in scope, would be useful to as broad a user base as possible. The goal was to design a software structure that was flexible for future needs while allowing for quick design comparison and  96  A.1. Introduction evaluation prior to any time-consuming fabrication work. The software system as designed is able to be implemented by a chip designer who is not necessarily trained in programming or mathematics to use it to help model and predict the performance of bead-based absorption or adsorption assays, in microfluidic chambers of varying geometry and design  A.1  Introduction  The motivation for developing a finite element analysis pipeline over an analytical approach aside from that already mentioned, can be illustrated by the following simple example. Suppose a bead with koff = 0 (perfect absorber) sits suspended in the middle of a well-mixed solution of concentration C(t) and volume V and that the rate of capture of molecules on the bead surface is diffusion-limited. The bead has a radius R and is surrounded by a constant width linear depletion region of diameter 2a. Such a system is illustrated schematically in Figure A.1. Depletion Region - Diameter 2a  Bead - Radius R  (a) Model diagram. A perfectly absorbing bead of radius R is surrounded by a depletion region of diameter 2a.  C(t1 ) = C1 C(t2 ) = C2 (b) Ideal concentration profile. The concentration outside the depletion region is constant and drops linearly within the depletion region to zero at the bead surface.  Figure A.1: Diagram for the analytical single bead model of a perfect absorber The solution for bulk concentration in this system (assuming a = 2R) describing C(t) at time t is shown in Equation A.1 where C0 is the bulk concentration at t = 0, D is the diffusion constant for the molecular species in question, and all other parameters are previously defined. The solution 97  A.1. Introduction is displayed graphically in Figure A.2. Full derivation details can be found in Appendix B.2. C(t) = C0 exp −  1.0  16πRD t V  (A.1)  Solution Concentration vs. Time  Normalized Concentration in Solution  0.8  0.6  0.4  0.2  0.0 0  5  10 Time (min)  15  20  Figure A.2: Concentration in a well-mixed solution over time with a single 5µm diameter bead acting as a perfect absorber in a 2nL chamber. Diffusion constant of 5 · 10−11 m2 /s. For this model a = 2R. 1 16πRD = gives τ ≈ 318s = 5.3min for τ V the well-mixed solution to decrease to 1/e of the initial concentration. This Defining a time constant  system describes the scenario of a cell instantaneously releasing a burst of molecules that are well-mixed to some uniform concentration C0 and then captured by the bead over time while remaining mixed. The derivation is trivial but the model has many non-physical assumptions. In the microfluidic devices designed, the capture occurs in a diffusion-limited regime with no forced convection (by design) meaning that the solution is never well-mixed. Imaging experiments have shown cell permeabilization to occur within a few seconds of exposure to mild lysis buffers and thus an instantaneous release may be an accurate enough approximation for an incubation-capture model; 98  A.1. Introduction however, since the system is not well-mixed, the concentration around the bead is not uniform at time zero, and in reality a single bead is not suspended within the center of the chamber. The analytical model may be improved slightly by extending it to incorporate a depletion region that increases in diameter as analyte capture progresses and the concentration is depleted. Setting up a similar system as before and solving analytically yields Equation A.2 (from [94]): A(r, t) = 1 −  1 1 − erf r  r−1 √ 2 t  (A.2)  In Equation A.2, A is the normalized concentration to initial concentration A0 , and r is the radial distance from (0,0) normalized to the bead radius, assuming a spherical bead whose geometric center is coincident with (0,0). η 2 2 The error function erf(η) is defined as erf(η) ≡ √ e−x dx [95]. The π 0 model no longer makes the incorrect assumption of a well-mixed solution; however, such a model still assumes a single bead in suspension with pulsed release being well mixed initially and adds the additional assumption of an infinite container in order to obtain an analytical solution. To more accurately approximate analyte capture on multiple beads, analytical approaches have been described for “capture elements” arranged in a spatial array with a single bead or porous element per unit cell [94, 96, 97]. For example, in [94], as capture progresses, the depletion regions overlap, effectively altering the boundary conditions of a unit cell. However, these systems must be solved numerically and still incorporate several limiting assumptions such as an infinitely extending periodic plane in R2 and can only handle periodic geometries. Furthermore, the systems derived as such as not user-friendly to quick manipulations for purposes of designing an assay or characterizing a system. To overcome these limitations and provide an accurate quantitative model for analyte capture on beads while enabling users unskilled in numerical methods or programming to make design decisions, a MATLAB-COMSOL finite element method analysis pipeline was designed.  99  A.1. Introduction  A.1.1  Finite Element Method and COMSOL Background  The modelling was performed using COMSOL’s coefficient form partial differential equation (PDE) physics mode where the governing system equation and boundary conditions are entered manually. The general governing PDE for this physics is defined internally to COMSOL [56] as: ea  ∂2u ∂u + da + ∇ · (−c∇u − αu + γ) + β · ∇u + au = f 2 ∂t ∂t  (A.3)  In Equation A.3, ea is the mass coefficient, da is the damping coefficient, c is the diffusion coefficient (a tensor that is isotropic, anisotropic, symmetrical, or diagonal), α is the conservative flux convection coefficient, γ is the conservative flux source, β is the convection coefficient, a is the absorption coefficient, and f is the source term. For the purpose of explaining the details of the construction of the model around this PDE, a 2-dimensional geometry is used here for illustrative purposes. An arbitrary polygon is constructed over which Equation A.3 is valid and is constrained to boundary and initial conditions [66]. This construction is shown in Figure A.3.  ΓN  ΓD Ω  Figure A.3: Arbitrary geometry in two dimensions for the coefficient form PDE with domain Ω in the plane R2 . The boundaries ΓN and ΓD represent Neumann and Dirichlet boundary conditions in thin and bold-face lines, respectively (Some sources use the notation of ∂ΩN and ∂ΩD ).  100  A.1. Introduction From Figure A.3, the problem can then be defined as follows:  ∂u ∂2u   + da + ∇ · (−c∇u − αu + γ) + β · ∇u + au = f, e a  2  ∂t ∂t    n · (c∇u + αu − γ) + qu = g − hT µ,   hu = r,      u(x, y) = fi (x, y),  in Ω on ΓN on ΓD at t = 0 (A.4)  The boundary conditions as written are taken from the notation used by COMSOL [57–59], where µ is a Lagrange multiplier, discussed in Section A.1.3, n is the unit normal vector, q is a boundary absorption term, g is a boundary source term, h is a scaling factor (not necessarily a scalar), and r is used to define Dirichlet boundary conditions. For a simple diffusion system with no convection and a damping coefficient of da = 1, Equation A.3 can be reduced to:  ∂u − ∇ · (c∇u) = 0 ∂t  (A.5)  More generally, this can be written as: ∂u(x, y, t) = ∇ · [D(u, x, y)∇u(x, y, t)] ∂t  (A.6)  Where u(x, y, t) is the concentration of species u at position (x, y) and time t in the constructed two-dimensional system, and D is the diffusion coefficient (as a function of concentration u, and position (x, y)). It is then further assumed that the diffusion coefficient D is constant and thus Equation A.6 reduces to the form: ∂u(x, y, t) = D∇2 u(x, y, t) ∂t  (A.7)  Next, it is assumed that −n·D∇u = g1 on ΓN and u = g0 on ΓD . Specifically, for the modelling performed, it is assumed that g0 = 0; that is, the bead acts as a perfect absorber and every molecule that comes into contact with the surface is instantly and irreversibly bound to the bead. As the bead surface contains many more binding sites than there are bound molecules — even  101  A.1. Introduction after extended periods of secretion — this approximation should hold for experimental time-frames. From this, the problem as described in (A.4) can be more precisely defined as:  ∂u(x, y, t)   = D∇2 u(x, y, t)   ∂t    u = g0   −D∂n u = g1      u(x, y) = fi (x, y),  in Ω on ΓD  (A.8)  on ΓN at t = 0  In case structure (A.8) g0 is a continuous function, g1 may be discontinuous, and −D∂n u = −c∇u · n where n is the normal vector. The versatility of the finite-element-method is based on re-writing the problem in a variational form called the “weak” form [66]. This is in contrast to the “strong” form, which is the problem as explicitly stated in (A.8). The strong form of a problem requires continuity on the dependent variables in differentiation up to the order of the highest appearing in the PDE [65]. The weak form by contrast, is an integral form with weaker constraints on continuity since the differentiation is split between a “weighting” or “test function” and the dependent variable. Typically, solving complex geometries in the strong form is much more difficult than in the weak form [65]. In addition, when discretizing the system, the strong form does not guarantee the same number of equations as unknowns, whereas this is guaranteed by the weak form of the problem [62]. Because of this, COMSOL converts all problem statements into the weak form internally when discretizing the system; this is done using Galerkin’s method of weighted residuals [58]. For more details on the method the interested reader may consult [61–64, 98]. A brief and general explanation of the weak form methodology as implemented by COMSOL is given below with the two-dimensional diffusion equation used for example. The governing system equation, as defined in case structure (A.8) is first integrated over the domain Ω while multiplied by an arbitrary weighting, or  102  A.1. Introduction test function, v. ∂u − D∇2 u dΩ = 0 ∂t  v Ω  (A.9)  Rearranging: v Ω  ∂u dΩ = D ∂t  v∇2 u dΩ  (A.10)  Ω  The gradient product rule is applied to the right hand side to give: ∂u dΩ = D ∂t  v Ω  ∇(v∇u) dΩ − D  ∇v∇u dΩ  Ω  (A.11)  Ω  The divergence theorem is then applied to the first term on the right hand side: v Ω  ∂u dΩ = D ∂t  (v∇u) · n ˆ dΓ − D  ∇v∇u dΩ  Γ  (A.12)  Ω  The boundary Γ is then split into Neumann and Dirichlet boundaries: v Ω  ∂u dΩ = D ∂t  (v∇u) · n ˆ dΓD + D  (v∇u) · n ˆ dΓN ΓN  ΓD  −D  ∇v∇u dΩ (A.13) Ω  Since v is arbitrary, it can be set to v = 0 on Dirichlet boundary conditions to eliminate the unknown flux terms. This simplifies the equation to: v Ω  ∂u dΩ = D ∂t  (v∇u) · n ˆ dΓN − D ΓN  ∇v∇u dΩ  (A.14)  Ω  Note that −D∇u · n ˆ is known and given in (A.8) as g1 on ΓN boundaries. Therefore: v Ω  ∂u dΩ = − ∂t  v g1 dΓN − D ΓN  ∇v∇u dΩ  (A.15)  Ω  A final rearrangement yields: v Ω  ∂u dΩ + ∂t  ∇v∇u dΩ = 0  v g1 dΓN + D ΓN  (A.16)  Ω  103  A.1. Introduction Equation A.16 is the weak form of the problem, subject to the same conditions as (A.8). To discretize Equation A.16, a mesh is created over the geometry dividing it into discrete elements over which the function is approximated by linear combinations of basis functions. For details see [57, 58, 62, 64]. It is thus assumed that u can be approximated as [99]: u(x, y, t)  U (x, y, t) =  ui (t)Ni (x, y)  (A.17)  i  Where ui is the value of the test function at the mesh nodes and Ni are basis functions. Using Galerkin’s method [62] the assumption is made that U = v and thus (A.16) becomes: R=  ui (t)Ni (x, y) Ω  i  i  ∂ui (t) Ni (x, y) dΩ ∂t ui (t)Ni (x, y)g1 dΓN  + ΓN  +D  i  ui (t)∇Ni (x, y) Ω  i  ui (t)∇Ni (x, y) dΩ (A.18) i  Because U is an approximation of the differential equation and not an exact solution, Equation A.18 will not sum to zero but rather some residual R. The finite element method works to minimize this residual, and thus obtain a more accurate representation of u by U . To minimize Equation A.18 the residual is differentiated with respect to u at each node i and set to zero. Note that this generates a linear system of equations with the same number of equations as unknowns. This system can then be solved to give the value of the dependent variable at each node and interpolated points via the basis functions, thus solving for U and approximating u.  A.1.2  Calculating Flux - 1D Heat Diffusion Example  If the instantaneous flux over a boundary at every time point is known, the number of molecules accumulated on a given bead at a given time can be easily determined. This directly correlates to fluorescent intensity of 104  A.1. Introduction the bead, should a label be introduced, and when summing over all beads, to the percentage of captured molecules. The elements of interest for flux calculation in this modelling are Dirichlet elements on the surface of the beads, where the surface concentration is assumed to be zero with the beads acting as perfect absorbers. In the standard FEM approach [64] Dirichlet nodes are eliminated from the system of equations to be solved; flux at these boundaries is then estimated after solving the system of equations by taking the derivative of the basis function on the boundary element at the point or surface of interest and integrating as required. This approach is inaccurate however, and leads to large errors sensitive to the discretization of the domain as will be shown. The discretization process of Section A.1.1 and the need to introduce the Lagrange Multipliers (See Section A.1.3) to deal with this inaccuracy is best illustrated through the use of an example. The following PDE system and matrix construction is adapted and modified from [64] with permission from Zimmerman, describing heat conduction in a 1D rod. The governing equation for heat conduction in one dimension with an energy source Q is described by: d dx  k  dT dx  +Q=0    k = 3.3 J/o C · m · s     Q = 10 J/m · s   T |x=0 = 1o C      q|x=1 = 1.25 J/m2 · s  (A.19)  (A.20)  Applying the discretization process as shown in Section A.1.1 with 3 evenly spaced mesh elements (giving 4 nodes) to Equation A.19 with boundary conditions and parameters described in (A.20) yields the linear system [K][Tn ] = −[qm ] + [Qm ] with stiffness matrix K, vector of unknowns Tn at each node, and flux and source contributions. The resulting system of equations is  105  A.1. Introduction shown in (A.21):   10 −10  0  0    −10 20 −10 0   0 −10 20 −10  0 0 −10 10    T1      q|x=0       T2   0     T  =  0  3   −1.25 T4      1.7         3.3   +   3.3     1.7  (A.21)  Applying the Dirichlet boundary condition T1 = T |x=0 = 1 and ignoring the first row, the matrices can be rearranged into a system of 3 equations and 3 unknowns   20 −10  0    T2      13.3         20 −10   T3  =  3.3   −10 0 −10 10 T4 0.45 This can be solved using standard Gaussian elimination to give:   T1      1         T2   1.705   =   T   2.08   3    T4 2.125 Comparing this to the analytical solution:   T1      1         T2   1.7095   =   T   2.089   3    T4 2.135 It is seen that even with 3 elements, the nodal values are already within 3% of the analytical solution. However, the flux at boundary nodes is not so accurately modelled. The flux at x = 0 is calculated by evaluating the derivative of the linear basis function from node 1 to node 2. This is simply: −k · ∆T /∆x = 7J/m2 · s  106  A.1. Introduction Comparing this to the analytical solution of 8.75J/m2 · s shows a 20% error in boundary flux. As this is the only Dirichlet boundary node, all energy added or removed from the system to maintain steady-state conditions must occur through this node and such a large difference is clearly unphysical and unacceptable. To calculate a more accurate flux value, the system can be rearranged to solve directly for q by applying the Dirichlet boundary condition T1 = 1:   10 −10  0  0    1      q|x=0       −10    20 −10 0     T2  =  0     0 −10 20 −10     T3   0 −1.25 0 0 −10 10 T4      1.7         3.3  +    3.3     1.7  Then rearranging to solve for the unknowns:        −1 −10  0  0      0    0 −10 20 −10   0 0 −10 10 0  20 −10  q|x=0 T2 T3 T4      −8.3         13.3  =    3.3     0.45  Reducing this system of equations yields:   q|x=0        T2 T3 T4      −8.75         1.705  =    2.08     2.125  Note that the nodal temperatures match the previous solution but the flux at x = 0 now matches the analytical solution. Since there is only one Dirichlet node in this example, all energy added to or removed from the system must pass through this node and therefore solving the system of equations for flux will always yield an exact answer. With more dimensions, the flux will be satisfied in a global average sense (with varying levels of accuracy on the Dirichlet boundary over each individual mesh element) because of the integral form of the PDE in the weak formulation. 107  A.1. Introduction However, this process of rearranging breaks the symmetry and sparsity of the stiffness matrix. The method of FEA matrix construction typically yields sparse tridiagonal banded matrices for the stiffness matrix [64]. Linear systems of this type can be solved by the Thomas Algorithm from which the solution can be obtained in O(n) operations instead of the O(n3 ) operations required for traditional Gaussian elimination [100, 101]. Breaking this structure thus greatly increases the computational complexity and required modelling time. Likely for this reason, among others, COMSOL does not allow matrix rearrangement in the manner shown above; thus solving using the traditional approach means flux calculations must be done after nodal values are determined leading to large inaccuracies and global imbalance of flux at steady state. To solve this accuracy problem, Lagrange multipliers were incorporated on Dirichlet boundaries (See Section A.1.3).  A.1.3  Lagrange Multipliers for the Finite Element Method  The method of Lagrange multipliers is a strategy for finding minima or maxima of a non-linear function of one or more variables and subject to equality constraints [95, 102, 103]. In general, for a two-variable function f (x, y) subject to g(x, y) = b the method of Lagrange multipliers seeks a solution: ∇Λ(x, y, λ) = ∇[f (x, y) + λ(g(x, y) − b)] = 0  (A.22)  Which is solved analytically. Or if an analytical solution is not available, numerical approaches attempt to minimize the residual: R(x, y, λ) = |∇Λ(x, y, λ)| = |∇[f (x, y) + λ(g(x, y) − b)]|  (A.23)  For multiple constraints gi (x1 , ..., xn ) = bi the Lagrange multipliers are implemented as [104]: m  λi (gi (x) − bi )  Λ(x, λ) = f (x) +  (A.24)  i=0  108  A.1. Introduction For a system of equations A · x = f with stiffness matrix A and nodal values xi generated by the FEM, and subject to Dirichlet boundary conditions B · x = h, the Lagrange multipliers λi can be incorporated by adding the Lagrange terms in Equation A.24 into Equation A.7 then proceeding with weak form construction and discretization to Equation A.18. This process yields the matrix system [105]: A BT  x  B  λ  0  =  f h  The Dirichlet conditions are now incorporated explicitly into the system of equations. The Lagrange multipliers then adjust to satisfy the equations now subject to these constraints. For the one dimensional heat problem in Section A.1.2 this yields:   0 1 0 0 0    T1    −10 20 −10 0 0 0 0 0   0 −10 20 −10 0 0 0 0    0 0 −10 10 0 0 0 0   1 0 0 0 0 0 0 0    0 0 0 0 0 0 0 0   0 0 0 0 0 0 0 0  0 0 0 0 0 0 0 0                       T2    3.3    T3   3.3       0.4  T4  =    λ1    1       λ2    0    λ3   0   λ4 0  10 −10  0      1.65    (A.25)  Note that this approach is similar to the rearrangement to solve directly for flux in Section but matrix symmetry and sparsity are maintained with the Lagrange multiplier approach — though the strict banded structure is lost. An inspection of (A.25) shows that the Lagrange multipliers have units of flux; depending on the system in question, they give the energy or mass transfer required to maintain the Dirichlet node or face at it’s prescribed value. Thus the flux given by the Lagrange multipliers will be globally conservative for steady state systems and over a continuous boundary will be accurate in an average sense across each Dirichlet surface due to the weighted residuals approach of the FEM. 109  A.1. Introduction The accuracy of this approach can be easily compared to the basis function derivative approach of Section A.1.2 since an analytical solution for this problem is known (q|x=0 = −8.75). The system is constructed in COMSOL and the number of elements varied to produce Table A.1. The accuracy of the Lagrange multiplier approach is evident as an exact solution is provided even with a single finite element. To approach a similar accuracy by computing derivatives even for this simple system requires 1000 elements. The computational cost for more complex systems to achieve comparable accuracy is immense. These results are shown graphically in Figure A.4. # of Elements 1 2 3 4 5 6 7 8 9 10 15 20 25 30 35 40 45 50 100 200 500 1000  −k · dT /dx -3.75 -6.25 -7.08 -7.50 -7.75 -7.91 -8.03 -8.12 -8.19 -8.25 -8.41 -8.49 -8.55 -8.58 -8.60 -8.62 -8.64 -8.64 -8.69 -8.72 -8.74 -8.75  λ -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75 -8.75  Table A.1: Comparison of basis function derivatives vs. Lagrange multiplier values to compute boundary flux (with different numbers of finite elements for the 1D heat example)  110  A.1. Introduction  −9  Calculated Flux vs. Number of Finite Elements (Truncated)  −8  2  Flux at x=0 mJ s  −7 −6 −5 −4 −3  Analytical Solution FEM Evaluated Derivative 0  50  100 Number of Elements  150  200  Figure A.4: Number of elements vs. flux accuracy for the traditional finite element method. To illustrate the accuracy difference between the derivative-based flux calculation and the Lagrange multipliers for a real modelled system and emphasize the effect on the global conservation of flux in steady state, a simple one cell - two bead model was constructed in COMSOL. Figure A.5 illustrates the geometry used for this comparison. A single cell was modelled as a simple sphere with radius rC = 5µm with two beads for capture acting as perfect absorbers and with a radius of rB = 2.5µm. For the derivative-based flux calculation the transport of diluted species physics module was used. The bead’s surface boundary was set at a constant concentration c = 0 while the exterior walls of the box were set to a no flux condition. For the Lagrange multiplier approach, the system was set-up manually in PDE form. The cell was set to have a total flux of 300molecules/s divided evenly over it’s surface in both cases. Figure A.6 shows the resulting flux calculations for the derivative-based approach. It is easily seen that as the system approaches steady state the flux difference between the cell and the beads remains quite large, which 111  A.1. Introduction  Figure A.5: Depiction of the geometry used for the basic diffusion and capture model. Drawing units are in µm. is physically impossible. Thus the capture efficiency and quantification of surface coverage will be inaccurate, as evidenced by Table A.1. Figure A.7 shows the resulting flux calculations for the Lagrange multiplier approach using the same system and parameters as implemented for Figure A.6 but incorporating the Lagrange multipliers to satisfy the Dirichlet boundary conditions. It is seen that as the system approaches steady state the flux onto the surface of the beads increases to match the number of molecules secreted from the cell ensuring global conservation of mass and maintaining a physically realizable system.  112  A.1. Introduction  Figure A.6: Cell and bead flux vs. time for the model depicted in Figure A.5 for the derivative-based approach to calculating boundary flux.  113  A.1. Introduction  Figure A.7: Cell and bead flux vs. time for the model depicted in Figure A.5 for the Lagrange multiplier approach to calculating boundary flux  114  A.2. Mean & Total Intensity Modelling  A.2  Mean & Total Intensity Modelling  At early time points, because of the effects mentioned, the mean intensity can sometimes provide a more accurate measurement than total intensity, but the trend will always go towards favouring total intensity measurements for long time periods. Figure A.9 illustrates this effect; the solid lines indicate the ‘true’ mean and the dashed lines indicate the measured mean, assuming a threshold value of 5000 molecules required before a signal is measurable. In contrast, Figure A.10 shows the much smaller perturbations in accuracy over time, trending quickly towards a uniform value for all numbers of beads and similar across distributions; note again however the greater accuracy of the far distribution. Figure A.8 shows a graph of the percent difference of the mean and total intensity for the various configurations of beads with the difference calculated between chambers with 60 beads and chambers with 140 beads.  115  A.2. Mean & Total Intensity Modelling  116 Figure A.8: Percent difference between the theoretically predicted measured mean and total intensity between a low-end number of beads (60) and a high-end number of beads (140) for different configurations in a 150 µm square chamber.  A.2. Mean & Total Intensity Modelling To properly interpret Figure A.8 it is necessary to determine what an effective incubation time is as in early time-points the errors can be comparable between mean and total intensity measurements. Figure 3.3 shows the percentage of molecules released by the cell that are captured by the beads over time for the various configurations. To achieve a high percentage capture, it is seen that around 120 min of capture time are needed; this is validated experimentally as well. Referring back to Figure A.8, in all cases, after sufficient incubation time, the total intensity measurement vastly outperforms the mean intensity measurement. Of special note is the extremely low percentage difference after 100 min produced by the far configuration.  117  A.2. Mean & Total Intensity Modelling  118 Figure A.9: Mean intensity measurement both theoretical and predicted actual measurement (based on a threshold of detection of 5000 molecules). The dashed lines indicate the predicted measurement value and the solid lines indicate the true value.  A.2. Mean & Total Intensity Modelling  119 Figure A.10: Total intensity measurement both theoretical and predicted actual measurement (based on a threshold of detection of 5000 molecules). The dashed lines indicate the predicted measurement value and the solid lines indicate the true value.  Appendix B  Derivations B.1  Concentration in a Box  The following derivation finds the solution to the system depicted pictorially in Figure B.1. The system in question is a sealed box, apart from one end, with a reagent mix flowing over the open side; it is desired to accurately determine the concentration profile in the box as a function of position and time. This derivation assumes that reagents are flowing over the open side of the box at such a rate that a molecule lost to diffusion within the box is instantly replenished; thus the concentration along this face is constant. It is further assumed that convection within the box is negligible and that all mixing happens via diffusion. Since the system is solved in one dimension, corner and edge effects are neglected. x=0  x=L  φ(0, t) = C0  ∂φ =0 ∂x  Figure B.1: System Definition  The system equation, as described by Fick’s second law is: ∂φ ∂2φ =D 2 ∂t ∂x  (B.1)  120  B.1. Concentration in a Box  In (B.1) φ is the concentration of a reagent with diffusion constant D and x and t are position and time variables, respectively. The boundary conditions, as depicted in Figure B.1 are as follows with C0 the concentration at x = 0: φ(0, t) = C0  (B.2)  ∂φ |x=L = 0 ∂x  (B.3)  The system is also subject to the initial condition φ(x, 0) = Ci , with Ci the initial concentration throughout the box (Restrained by Ci ≤ C0 ): φ(x, 0) = Ci  (B.4)  To begin solving this system, it is assumed that φ can be written in the following form: φ(x, t) = W (x, t) + V (x)  (B.5)  At steady state (t → ∞) W (x, t) → 0 since concentration cannot increase indefinitely. Therefore: φ(x, t)|t→∞ = V (x) = C0  (B.6)  ∂2φ ∂φ |t→∞ = 0 = ⇒ V (x) = Ax + B ∂t ∂x2  (B.7)  Applying (B.2) to (B.7):  121  B.1. Concentration in a Box V (0) = A(0) + B = C0  (B.8)  V (0) = B = C0  (B.9)  ∴ V (x) = Ax + C0  (B.10)  Applying (B.3) to (B.7): ∂φ ∂V |t→∞ = =A ∂x ∂x ∂V |x=L = 0 = A ∂x ∴ V (x) = C0  (B.12)  ∴ φ(x, t) = W (x, t) + C0  (B.14)  (B.11)  (B.13)  Rearranging this solution for W (x, t) yields: W (x, t) = φ(x, t) − C0  (B.15)  This rearrangement has new boundary conditions. The first can be derived by application of (B.2): W (0, t) = φ(0, t) − C0 = C0 − C0 = 0  (B.16)  W (0, t) = 0  (B.17)  The second can be derived using (B.3): ∂ [W (x, t) = φ(x, t) − C0 ] ∂x ∂W ∂φ = ∂x ∂x ∂φ ∂W |x=L = 0 = |x=L ∂x ∂x  (B.18) (B.19) (B.20)  122  B.1. Concentration in a Box  The new boundary conditions for W (x, t) are thus: W (0, t) = 0  (B.21)  ∂W |x=L = 0 ∂x  (B.22)  Next, it is assumed that W (x, t) is separable in x and t: W (x, t) = X(x)T (t)  (B.23)  The original differential equation (B.1) can now be re-defined in terms of W . Since V is linear in x, this reduces to: ∂2W ∂W =D ∂t ∂x2  (B.24)  Simplifying (B.24), rearranging, and introducing the constant eigenvalue λ yields: T X = DT X  (B.25)  T X = =λ DT X  (B.26)  There are then three possible cases for λ: λ=0 λ = σ2 λ = −σ 2 Case 1 (λ = 0): Starting with X: 123  B.1. Concentration in a Box X = λX  (B.27)  X = 0 ⇒ X(x) = Ax + B  (B.28)  New conditions are needed to solve this system. From (B.21): W (0, t) = X(0)T (t) = 0  (B.29)  For non-trivial solutions, X(0) = 0. From (B.22): ∂W |x=L = X (L)T (t) = 0 ∂x  (B.30)  For non-trivial solutions, X (L) = 0. Therefore: X(x) = Ax + B  (B.31)  X(0) = B = 0  (B.32)  X (L) = A = 0  (B.33)  Therefore X(x) = 0. The solution is trivial and therefore λ = 0. Case 2 (λ = σ 2 ): Again, starting with X: X − σ2X = 0  (B.34)  The characteristic equation for such a system is r2 − σ 2 = 0 with roots of ±σ. ∴ X(x) = c1 eσx + c2 e−σx  (B.35)  124  B.1. Concentration in a Box  New conditions are needed to solve this system. From (B.21): W (0, t) = X(0)T (t) = 0  (B.36)  For a non-trivial solution X(0) = 0. Now, solving for X(x) using (B.22): X(0) = c1 + c2 = 0 ⇒ c1 = −c2 σx  ∴ X (x) = c1 σ(e  σL  X (L) = c1 σ(e  −σx  +e  +e  )=0  −σL  ) = 0 ⇒ c1 = 0 = c2  (B.37) (B.38) (B.39)  Again, this is a trivial solution so λ = −σ 2 . Case 3 (λ = −σ 2 ): Starting with X again: X + σ2X = 0  (B.40)  The characteristic equation for such a system is r2 + σ 2 = 0 with roots of √ ±iσ, i = −1. ∴ X(x) = c1 cos(σx) + c2 sin(σx)  (B.41)  From (B.21) and (B.22): X(0) = c1 = 0 (2n − 1)π , n = 1, 2, 3, ... 2L (2n − 1)πx ∴ X(x) = c2 sin , n = 1, 2, 3, ... 2L  X (x)|x=L = c2 σ cos(σL) = 0 ⇒ σ =  (B.42) (B.43) (B.44)  125  B.1. Concentration in a Box  Now, looking at T: T + σ2T = 0 2  (B.45)  T = −σ DT ⇒ T = c1 e  −σ 2 Dt  (B.46)  Combining (B.44) and (B.46) by (B.23) where σ is defined by (B.43) and letting c2 c1 = bn : W (x, t) = X(x)T (t)  (B.47)  (2n − 1)πx exp − 2L  W (x, t) = bn sin  (2n − 1)πx 2L  2  Dt , n = 1, 2, 3, ... (B.48)  Each value of n yields a linearly independent solution for W (x, t). Therefore, the full solution for W (x, t) is a linear combination of each solution from n = 1 → ∞. This yields: ∞  W (x, t) =  bn sin n=1  (2n − 1)πx exp − 2L  (2n − 1)πx 2L  2  Dt  (B.49)  Substituting (B.49) into (B.5) and applying (B.6) yields: ∞  bn sin  φ(x, t) = C0 + n=1  (2n − 1)πx exp − 2L  (2n − 1)πx 2L  2  Dt  (B.50)  bn can be found by applying the initial condition (B.4):  126  B.2. Analytical Model for Perfect Absorber Bead Capture ∞  φ(x, 0) = C0 +  bn sin  (2n − 1)πx = Ci 2L  (B.51)  bn sin  (2n − 1)πx = Ci − C0 2L  (B.52)  n=1 ∞ n=1  Using the Fourier Inversion Theorem: L 2 (2n − 1)πx dx (Ci − C0 ) sin L 2L 0 2 (2n − 1)π 2 2L bn = − (Ci − C0 ) cos + (Ci − C0 ) L 2 L (2n − 1)π  bn =  (B.53) (B.54)  Simplifying: 2 (2n − 1)π − (Ci − C0 ) cos =0∀n L 2 4(Ci − C0 ) ∴ bn = π(2n − 1)  (B.55) (B.56)  Substituting (B.56) into (B.50) yields the full solution: ∞  φ(x, t) = C0 + n=1  4(Ci − C0 ) (2n − 1)πx −(2n − 1)2 π 2 Dt sin exp π(2n − 1) 2L 4L2 (B.57)  B.2  Analytical Model for Perfect Absorber Bead Capture  The governing equation for such a system as depicted in Figure B.2 is Fick’s first law of diffusion in polar coordinates: JD = −D∇C  (B.58)  127  B.2. Analytical Model for Perfect Absorber Bead Capture  Depletion Region - Diameter 2a C(t1 ) = C1 C(t2 ) = C2  Bead - Radius R  (a) Model diagram. A perfectly absorbing bead of radius R is surrounded by a depletion region of diameter 2a.  (b) Ideal concentration profile. The concentration outside the depletion region is constant and drops linearly within the depletion region to zero at the bead surface.  Figure B.2: Diagram for the analytical single bead model of a perfect absorber The initial condition for this system is: C(t = 0) = C0  (B.59)  Based on the symmetry of a spherical system, (B.58) can be reduced to: JD = −D  ∂C ∂r  (B.60)  To find the total molecular flux across the depletion region surface, (B.60) is integrated over the surface area. I=  JD · dA  Using the definition of the first derivative in polar coordinates  (B.61)  ∆C ∆r  the  equation can be simplified as follows:  128  B.2. Analytical Model for Perfect Absorber Bead Capture  I = JD  dA  (B.62)  −D(C(t) − 0) 4πa2 (a − R) 4πa2 DC(t) I=− (a − R) I=  (B.63) (B.64)  Integrating the total flux over time and subtracting this from the initial concentration in bulk solution, yields the concentration in bulk solution at time t: t  C(t) = C0 − 0 t  C(t) = C0 − 0  |I| dt V 4πa2 DC(t) dt V (a − R)  (B.65) (B.66)  Differentiating both sides of (B.66) with respect to t yields: 4πa2 DC(t) d[C(t)] =− dt (a − R)  (B.67)  Assume that the solution to C(t) takes the form: C(t) = c1 e−c2 t  (B.68)  Constant c1 can be solved for by applying boundary condition (B.59) to (B.68): C(0) = C0 = c1  (B.69)  Substituting (B.68) into (B.67) with c1 = C0 yields (B.70), which can be reduced to solve for the constant c2 :  129  B.2. Analytical Model for Perfect Absorber Bead Capture 4πa2 DC0 e−c2 t V (a − R) 2 4πa D c2 = V (a − R) 4πa2 D ∴ C(t) = C0 exp − t V (a − R) −C0 c2 e−c2 t = −  (B.70) (B.71) (B.72)  If a = 2R then (B.72) reduces to (B.73): C(t) = C0 exp −  16πRD t V  (B.73)  Solution (B.73) is depicted graphically in Figure B.3.  1.0  Solution Concentration vs. Time  Normalized Concentration in Solution  0.8  0.6  0.4  0.2  0.0 0  5  10 Time (min)  15  20  Figure B.3: Concentration in a well-mixed solution over time with a single 5µm bead acting as a perfect absorber in a 2nL chamber. Diffusion constant of 5 · 10−11 m2 /s. For this model a = 2R.  130  Appendix C  Materials & Methods C.1  LabView Interfaces  Figure C.1: LabView interface used for timing and solenoid/valve control.  131  C.1. LabView Interfaces  132 Figure C.2: Main LabView interface used for microscope control and imaging.  Appendix D  Coding D.1  1  2  3 4 5 6 7  Bead and Cell Chamber Configuration  % Script for creating a 'position file' with x−y coordinates of... beads and % a cell within a microfluidic chamber distributed according to... various % chosen parameters % % % Date: February 2013 % Author: Daniel J. Lewis  8 9  10  11 12 13 14  % Close all open figures, clear workspace, and clear command ... window %... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−... close all clear all clc %... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−...  15 16 17  18 19  % Set parameters %... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−... r = 4.90/2; % bead radius in micron rshift = 0.03; % slight increase in radius so beads don't touch  133  D.1. Bead and Cell Chamber Configuration 20 21 22 23 24 25 26 27  28  29  30 31 32  r = r+rshift; cellr = 15/2; % cell radius in micron xmax = 150; % chamber width in micron ymax = 150; % chamber length in micron N = 100; % total number of beads in the chamber s = 100; % max number of beads on the 'low occupancy' side smin = 70; % min number of beads on the 'low occupancy' side n = randi([smin s]); % create a random number of beads on the '... low' side low occ frac = 0.5; % ratio of low occ. length to total chamber... length % low occ len is 'length' of low occ. side, relative to far ... wall low occ len = xmax* low occ frac; xdiv = xmax − low occ len; %... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−...  33 34  %... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−...  35  cellpos = [0,0]; % initialize cell position vector cellpos(1) = xmax/2; % set the x−position of the cell cellpos(2) = 0+cellr+1; % set the y−position of the cell, off ... from wall cellposrange x = [0,0]; % initialize vector for cell occupied ... space %{ cellposrange x and cellposrange y set the x and y end−points ... within which the cell occupies space in the x−y plane and therefore ... where other objects are not allowed to reside. It is used to remove un−... physical locations from the possible bead location array. %} cellposrange x(1) = cellpos(1)−cellr−r; cellposrange x(2) = cellpos(1)+cellr+r; cellposrange y(1) = 0;  36 37  38  39 40  41  42  43 44 45 46 47  134  D.1. Bead and Cell Chamber Configuration 48  cellposrange y(2) = cellpos(2)+cellr+r;  49 50 51 52 53 54 55 56 57  x = r:2*r:xmax−r; shiftx = (xmax−max(x)−r)/2; x = x+shiftx; y = r:2*r:ymax−r; shifty = (ymax−max(y)−r)/2; y = y+shifty; [X,Y] = meshgrid(x,y); rs = find(x>xdiv,1,'first');  58 59 60  remove x val = remove y val =  find(cellposrange x(1)<x&x<cellposrange x(2)); find(cellposrange y(1)<x&x<cellposrange y(2));  61 62 63 64 65 66 67 68 69  remove x = zeros(length(remove x val),1); remove y = zeros(length(remove y val),1); for i = 1:length(remove x val); remove x(i) = x(remove x val(i)); end for i = 1:length(remove y val); remove y(i) = y(remove y val(i)); end  70 71 72  73  74  75 76 77 78 79 80 81 82 83 84  %{ The for−loop structure below steps through all the combinations... of X−Y coordinates and tests them against the un−allowed positions ... within the cell area. If the X−Y point lies within this area the pos−... array holds a value of [0,0] else it holds [X(i,j),Y(i,j]) %} pos = cell(length(x),length(y)); remhold = [0,0]; % initialize holding vector for X Y data for i = 1:length(x) for j = 1:length(y) remhold(1) = X(i,j); % hold x−val for current iteration remhold(2) = Y(i,j); % hold y−val for current iteration if max(remhold(1)==remove x)==1 if max(remhold(2)==remove y)==1  135  D.1. Bead and Cell Chamber Configuration pos{i,j} = [0,0]; else pos{i,j} = [X(i,j),Y(i,j)]; end  85 86 87 88  else  89  pos{i,j} = [X(i,j),Y(i,j)];  90  end  91  end  92 93  end  94 95  96  97 98 99 100 101 102 103 104 105 106  % The loop below steps through the pos cell array to create a ... matrix. You % end up with an Nx2 matrix composed of the columns of pos ... concatenated % together sizepos = size(pos); k = 1; posM = []; for i = 1:sizepos(1) for j = 1:sizepos(2) posM = cat(1,posM,pos{j,i}); k = k+1; end end  107 108 109  110  111  112  113  114 115 116 117  %{ This first little for−loop structures just creates a row−... position vector (a2) with duplicates from a removed (since both the x and y ... columns of posM that are set to 0 will show up in a. b gives the column ... number and a gives the row number with all of column 1 coming before ... column 2 so this loop just truncates a before the positions start ... duplicating %} [a,b] = find(posM==0); bfirs = find(b==1); a2 = zeros(length(bfirs),1);  136  D.1. Bead and Cell Chamber Configuration 118 119 120  for i = 1:length(bfirs) a2(i)=a(i); end  121 122  123  124  125  126  127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143  % a2 is set to integer type so that integer division can be ... performed. The % for loop then checks to see which column within the pos cell ... array the % row is equivalent too (since the columns were all ... concatenated for posM) % by doing a division check to see how many column equivalents ... (based on % number of rows) have been 'passed'. a2col is a vector ... containing the % column number for each row in a2 (a). a = a2; a2 = int32(a2); for i = 1:length(a2); if a2(i)>30 quo = idivide(a2(i),int32(sizepos(1))); remain = rem(a2(i),int32(sizepos(1))); if remain>0 a2(i)=a2(i)−quo*sizepos(1); a2col(i) = quo+1; end if remain==0 a2(i)=sizepos(1); a2col(i) = quo; end end end  144 145  146  147  % a2 and a2col are concatenated to form a row/column pair for ... the pos cell % array entry to be removed (since it lies within the cell area... ) a2 = cat(2,a2,a2col');  148 149 150  [row,col] = size(pos); row = 1:row;  137  D.1. Bead and Cell Chamber Configuration 151  col = 1:col;  152 153 154 155 156  maxl = max(row)*max(col); k = 1; kp = 1; j = 1;  157 158  holdv = zeros(maxl,2);  159 160  161  162  163 164 165 166 167 168 169 170 171 172 173 174  % this for−loop creates a matrix (hold−v) whose first column ... steps through % row numbers of pos (1,2,3..# rows) and repeats for the number... of columns. % The second column steps through column numbers every time the... max number % of rows is reached(1,1,1.. for # rows.. 2,2,2..) for i = 1:maxl; holdv(i,1) = row(j); holdv(i,2) = col(k); kp = kp+1; j = j+1; if kp==length(row)+1 k = k+1; kp = 1; j = 1; end end  175 176 177  178 179 180  181  holdv mem = holdv; % the rows defined by a above that refer to positions that lie ... within the % cell are then removed holdv = removerows(holdv,a); % the rows are then sorted by the first column in ascending ... order holdv = sortrows(holdv);  182 183  184  % mxl finds the last row that resides within the 'high ... occupancy' side of % the chamber  138  D.1. Bead and Cell Chamber Configuration 185  mxl = find(holdv(:,1)<=rs,1,'last');  186 187  188  189  190 191 192 193 194 195 196 197  198 199 200 201 202 203 204  % loc is a random sample of integers from 1 to mxl with N−n ... element (the % number of cells on the 'high occupancy' side. matloc is a ... new matrix % which takes the cell array position data from holdv ... corresponding to the % randomly selected rows loc = randsample(mxl,N−n); matloc = zeros(N−n,2); for i = 1:(N−n) matloc(i,1) = holdv(loc(i),1); matloc(i,2) = holdv(loc(i),2); end % loc2 and matloc2 does the same thing as loc but for the 'low ... occupancy' % side loc2 = randsample(mxl+1:length(holdv),n); matloc2 = zeros(n,2); for i = 1:n matloc2(i,1) = holdv(loc2(i),1); matloc2(i,2) = holdv(loc2(i),2); end  205 206  207 208  % the two location matrices are concatenated together to ... address the pos % cell array and get the final X−Y coordinate data matloc = cat(1,matloc,matloc2);  209 210  211  212 213 214 215 216 217  % this for loop uses the row−column identifiers from matloc and... pulls the % X−Y position data from the respective position in the pos ... cell array posloc = zeros(N,2); for i = 1:N poshold = pos{matloc(i,1),matloc(i,2)}; posloc(i,1) = poshold(1); posloc(i,2) = poshold(2); end  139  D.1. Bead and Cell Chamber Configuration 218 219  220  221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242  % the lines and nested if else statement below are not ... currently being % used; they separate out the x−y position data based on which ... 'column' of % y they belong to. sortpos = sortrows(posloc); h = 1; j = 1; posbycol = cell(1,1); for i = 1:length(sortpos) hold j(i) = j; if i==1 posbycol{h,j} = sortpos(i,:); h = h+1; else if sortpos(i,1)==sortpos(i−1,1) posbycol{h,j} = sortpos(i,:); h = h+1; else h = 1; j = j+1; posbycol{h,j} = sortpos(i,:); end end end hold j = hold j';  243 244 245 246 247 248 249  fig1 = scatter(posloc(:,1),posloc(:,2),32*2,'k','filled'); hold on; scatter(cellpos(1),cellpos(2),64*8,'r','filled'); grid off; xlim([0 xmax]); ylim([0 ymax]);  250 251  r = r−rshift;  252 253 254  save bead config.out r cellr xmax ymax N cellpos posloc −ASCII saveas(gcf,'beadfig','png');  140  D.2. Bead Capture Model Construction  D.2  1  Bead Capture Model Construction  function [t1] = beadmodel(num1,num2,num3)  2 3  % Model initially exported on Feb 7 2013, 16:38 by COMSOL 4... .0.0.982.  4 5  tic;  6 7 8  import com.comsol.model.* import com.comsol.model.util.*  9 10  model = ModelUtil.create('Model');  11 12  model.modelPath(pwd);  13 14  model.modelNode.create('mod1');  15 16  model.geom.create('geom1', 3);  17 18  model.mesh.create('mesh1', 'geom1');  19 20  model.physics.create('c', 'CoefficientFormPDE', 'geom1', {'u'})... ;  21 22 23  model.study.create('std1'); model.study('std1').feature.create('time', 'Transient');  24 25  model.geom('geom1').lengthUnit([native2unicode(hex2dec('00b5'),... 'Cp1252') 'm']);  26 27 28  readpath = strcat('bead config',num1,' ',num2,' ',num3,'.out'); dat in = dlmread(readpath);  29 30 31 32 33 34  r = dat in(1,1); cellr = dat in(2,1); xmax = dat in(3,1); ymax = dat in(4,1); N = dat in(5,1);  141  D.2. Bead Capture Model Construction 35 36 37 38 39  cellpos(1) = dat in(6,1); cellpos(2) = dat in(6,2); posloc = dat in(7:length(dat in),:); area per bead = 4*pi()*rˆ2; % area per bead in micron squared total bead area = area per bead *N; % total exposed bead area;  40 41 42  43  44  45  model.geom('geom1').feature.create('blk1', 'Block'); model.geom('geom1').feature('blk1').setIndex('size', num2str(... xmax), 0); model.geom('geom1').feature('blk1').setIndex('size', num2str(... ymax), 1); model.geom('geom1').feature('blk1').setIndex('size', num2str(... xmax), 2); model.geom('geom1').run('blk1');  46 47  model.view('view1').set('transparency', 'on');  48 49 50 51  52  53  54  model.geom('geom1').feature.create('sph1', 'Sphere'); model.geom('geom1').feature('sph1').set('r', num2str(cellr)); model.geom('geom1').feature('sph1').setIndex('pos', num2str(... cellpos(1)), 0); model.geom('geom1').feature('sph1').setIndex('pos', num2str(... cellpos(2)), 1); model.geom('geom1').feature('sph1').setIndex('pos', num2str(... cellr), 2); model.geom('geom1').run('sph1');  55 56 57 58 59 60 61 62  63  64  65  for i = 1:length(posloc) fnum = i+1; fnum = num2str(fnum); featname = strcat('sph',fnum); model.geom('geom1').feature.create(featname,'Sphere'); model.geom('geom1').feature(featname).set('r',num2str(r)); model.geom('geom1').feature(featname).setIndex('pos',... num2str(posloc(i,1)),0); model.geom('geom1').feature(featname).setIndex('pos',... num2str(posloc(i,2)),1); model.geom('geom1').feature(featname).setIndex('pos',... num2str(r),2); model.geom('geom1').run(featname);  142  D.2. Bead Capture Model Construction holdi = i;  66 67  end  68 69 70  model.geom('geom1').run; model.geom('geom1').runAll;  71 72 73 74  disp('Build Geometry: Pass') toc; tic;  75 76 77 78  oNames = model.geom('geom1').objectNames(); upDown = model.geom('geom1').getUpDown; upDown = int32(upDown');  79 80  domain = cell(1,length(upDown));  81 82 83 84  for i = 1:length(domain) domain{upDown(i,2)} = cat(2,domain{upDown(i,2)},i); end  85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105  sortpos = sortrows(posloc); xlessrows = find((sortpos(:,1)−r)<(cellpos(1)−cellr)); selcount = 3; bndset beads = []; for i = 2:(max(xlessrows)+1); bndset beads = cat(2,bndset beads,domain{i}); selnum = i+1; selnum = strcat('sel',num2str(selnum)); model.selection.create(selnum); model.selection(selnum).geom(2); model.selection(selnum).set(domain{i}); selcount = selcount+1; end for i = (max(xlessrows)+3):length(domain); bndset beads = cat(2,bndset beads,domain{i}); selnum = selcount; selnum = strcat('sel',num2str(selnum)); model.selection.create(selnum); model.selection(selnum).geom(2); model.selection(selnum).set(domain{i});  143  D.2. Bead Capture Model Construction 106 107 108  selcount = selcount+1; end bndset cell = domain{max(xlessrows)+2};  109 110 111 112  disp('Create Selection: Pass') toc; tic;  113 114 115 116  117  118  119  120  121  122  123  124  125  126  127  128  model.material.create('mat1'); model.material('mat1').name('Water'); model.material('mat1').materialModel('def').set('... dynamicviscosity', 'eta(T[1/K])[Pa*s]'); model.material('mat1').materialModel('def').set('... ratioofspecificheat', '1.0'); model.material('mat1').materialModel('def').set('... electricconductivity', '5.5e−6[S/m]'); model.material('mat1').materialModel('def').set('heatcapacity',... 'Cp(T[1/K])[J/(kg*K)]'); model.material('mat1').materialModel('def').set('density', 'rho... (T[1/K])[kg/mˆ3]'); model.material('mat1').materialModel('def').set('... thermalconductivity', 'k(T[1/K])[W/(m*K)]'); model.material('mat1').materialModel('def').set('soundspeed', '... cs(T[1/K])[m/s]'); model.material('mat1').materialModel('def').func.create('eta', ... 'Piecewise'); model.material('mat1').materialModel('def').func('eta').set('... funcname', 'eta'); model.material('mat1').materialModel('def').func('eta').set('... arg', 'T'); model.material('mat1').materialModel('def').func('eta').set('... extrap', 'constant'); model.material('mat1').materialModel('def').func('eta').set('... pieces', {'273.15' '413.15' '1.3799566804−0.021224019151*T... ˆ1+1.3604562827E−4*Tˆ2−4.6454090319E−7*Tˆ3+8.9042735735E... −10*Tˆ4−9.0790692686E−13*Tˆ5+3.8457331488E−16*Tˆ6'; '413.15... ' '553.75' '0.00401235783−2.10746715E−5*Tˆ1+3.85772275E−8*T... ˆ2−2.39730284E−11*Tˆ3'}); model.material('mat1').materialModel('def').func.create('Cp', '... Piecewise');  144  D.2. Bead Capture Model Construction 129  130  131  132  133  134  135  136  137  138  139  140  141  142  143  144  145  146  model.material('mat1').materialModel('def').func('Cp').set('... funcname', 'Cp'); model.material('mat1').materialModel('def').func('Cp').set('arg... ', 'T'); model.material('mat1').materialModel('def').func('Cp').set('... extrap', 'constant'); model.material('mat1').materialModel('def').func('Cp').set('... pieces', {'273.15' '553.75' '12010.1471−80.4072879*Tˆ1+0... .309866854*Tˆ2−5.38186884E−4*Tˆ3+3.62536437E−7*Tˆ4'}); model.material('mat1').materialModel('def').func.create('rho', ... 'Piecewise'); model.material('mat1').materialModel('def').func('rho').set('... funcname', 'rho'); model.material('mat1').materialModel('def').func('rho').set('... arg', 'T'); model.material('mat1').materialModel('def').func('rho').set('... extrap', 'constant'); model.material('mat1').materialModel('def').func('rho').set('... pieces', {'273.15' '553.75' '838.466135+1.40050603*Tˆ1−0... .0030112376*Tˆ2+3.71822313E−7*Tˆ3'}); model.material('mat1').materialModel('def').func.create('k', '... Piecewise'); model.material('mat1').materialModel('def').func('k').set('... funcname', 'k'); model.material('mat1').materialModel('def').func('k').set('arg'... , 'T'); model.material('mat1').materialModel('def').func('k').set('... extrap', 'constant'); model.material('mat1').materialModel('def').func('k').set('... pieces', {'273.15' '553.75' '−0.869083936+0.00894880345*T... ˆ1−1.58366345E−5*Tˆ2+7.97543259E−9*Tˆ3'}); model.material('mat1').materialModel('def').func.create('cs', '... Interpolation'); model.material('mat1').materialModel('def').func('cs').set('... funcname', 'cs'); model.material('mat1').materialModel('def').func('cs').set('... interp', 'piecewisecubic'); model.material('mat1').materialModel('def').func('cs').set('... extrap', 'const');  145  D.2. Bead Capture Model Construction 147  148  model.material('mat1').materialModel('def').func('cs').set('... table', {'273' '1403'; '278' '1427'; '283' '1447'; '293' '... 1481'; '303' '1507'; '313' '1526'; '323' '1541'; '333' '... 1552'; '343' '1555'; '353' '1555'; '363' '1550'; '373' '... 1543'}); model.material('mat1').materialModel('def').addInput('... temperature');  149 150 151 152 153 154 155  model.selection.create('sel1'); model.selection('sel1').geom(2); model.selection('sel1').set(bndset cell); model.selection.create('sel2'); model.selection('sel2').geom(2); model.selection('sel2').set(bndset beads);  156 157 158 159 160  model.variable.create('var1'); model.variable('var1').model('mod1'); model.variable('var1').set('secretion', '162.5[1/s]'); model.variable('var1').set('avogadro', '6.02e23[1/mol]');  161 162 163 164  model.cpl.create('intop1', 'Integration', 'geom1'); model.cpl('intop1').selection.geom('geom1', 2); model.cpl('intop1').selection.named('sel1');  165 166  model.physics('c').feature('cfeq1').set('c', 1, {'1' '0' '0' '0... ' '1' '0' '0' '0' '1'});  167 168  model.variable('var1').set('diffusion', '4.4e−11[mˆ2/s]');  169 170  171 172 173 174  175 176  177  model.physics('c').feature('cfeq1').set('c', 1, {'diffusion' '0... ' '0' '0' 'diffusion' '0' '0' '0' 'diffusion'}); model.physics('c').feature('cfeq1').set('f', 1, '0'); model.physics('c').feature.create('flux1', 'FluxBoundary', 2); model.physics('c').feature('flux1').selection.named('sel1'); model.physics('c').feature('flux1').set('g', 1, 'secretion/... intop1(1)/avogadro'); model.physics('c').feature.create('wc1', 'WeakConstraint', 2); model.physics('c').feature('wc1').set('constraintExpression', ... 1, 'u'); model.physics('c').feature('wc1').selection.named('sel2');  146  D.2. Bead Capture Model Construction 178  model.physics('c').feature('wc1').set('fieldVariableName', 1, '... flux');  179 180 181  182  183  184  185  186  187  188  189  190  191  192  193  model.mesh('mesh1').feature.create('ftet1', 'FreeTet'); model.mesh('mesh1').feature('ftet1').feature.create('size1', '... Size'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... custom', 'on'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hmaxactive', 'on'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hminactive', 'on'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hgradactive', 'on'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hcurveactive', 'on'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hnarrowactive', 'on'); model.mesh('mesh1').feature('ftet1').feature('size1').set('hmax... ', '15'); model.mesh('mesh1').feature('ftet1').feature('size1').set('hmin... ', '0.1'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hgrad', '2'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hcurve', '0.5'); model.mesh('mesh1').feature('ftet1').feature('size1').set('... hnarrow', '1'); model.mesh('mesh1').run;  194 195 196 197  disp('Build Mesh: Pass') toc; tic;  198 199 200 201 202  203  tmin = 0; tmax = 10000; tstep = 10; trange = strcat('range(',num2str(tmin),',',num2str(tstep),',',... ... num2str(tmax),')');  147  D.2. Bead Capture Model Construction 204  model.study('std1').feature('time').set('tlist', trange);  205 206 207 208 209 210 211 212  213 214 215 216 217 218  219 220  model.sol.create('sol1'); model.sol('sol1').feature.create('st1', 'StudyStep'); model.sol('sol1').feature('st1').set('study', 'std1'); model.sol('sol1').feature('st1').set('studystep', 'time'); model.sol('sol1').feature.create('v1', 'Variables'); model.sol('sol1').feature.create('t1', 'Time'); model.sol('sol1').feature('t1').set('tlist', 'range(0,25,3000)'... ); model.sol('sol1').feature('t1').set('plot', 'off'); model.sol('sol1').feature('t1').set('plotfreq', 'tout'); model.sol('sol1').feature('t1').set('probesel', 'all'); model.sol('sol1').feature('t1').set('probes', {}); model.sol('sol1').feature('t1').set('probefreq', 'tsteps'); model.sol('sol1').feature('t1').feature.create('fc1', '... FullyCoupled'); model.sol('sol1').feature('t1').feature.remove('fcDef'); model.sol('sol1').attach('std1');  221 222  model.sol('sol1').runAll;  223 224 225 226  disp('System Solved: Pass') toc; tic;  227 228 229 230 231 232 233 234 235 236 237 238  239 240  flux array = cell(1,1); for i = 1:(N+2) intnum = i; intnum = strcat('int',num2str(intnum)); selnum = i; selnum = strcat('sel',num2str(selnum)); model.result.numerical.create(intnum,'IntSurface'); model.result.numerical(intnum).selection.named(selnum); if i==1; model.result.numerical(intnum).set('expr', 'c.g u'); flux array{1,i} = model.result.numerical(intnum)... .getReal(); else model.result.numerical(intnum).set('expr','flux');  148  D.2. Bead Capture Model Construction 241  242 243 244  flux array{1,i} = model.result.numerical(intnum)... .getReal(); end end x time = linspace(tmin,tmax,length(flux array{1,1}));  245 246 247  248 249 250 251 252 253 254 255 256 257 258  259 260 261 262 263 264 265 266  267 268 269  % Calculating total and mean flux and captured molecules % ... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−... flux array quad = zeros(length(flux array{1,1})−1,N+2); flux array sum = zeros(length(flux array{1,1})−1,N+2); for i=1:(N+2) for j = 1:(length(flux array{1,1})−1) flux array Mhold = flux array{i}; flux array quad(j,i) = (x time(j+1)−x time(j))*... ((flux array Mhold(j+1)+flux array Mhold(j))/2); flux array sum(j,i) = sum(flux array quad(:,i)); end end flux array sum = flux array sum.*6.02e23; % total captured ... molecules tot perc capt = zeros(length(flux array{1,1})−1,1); for i = 1:(length(flux array{1,1})−1) tot perc capt(i) = flux array sum(i,2)/flux array sum(i,1); end tot intensity = zeros(length(flux array{1,1})−1,1); for i = 1:(length(flux array{1,1})−1) for j = 3:(N+2) tot intensity(i) = tot intensity(i)+flux array sum(i,j)... ; end end mean intensity = tot intensity./total bead area;  270 271 272 273 274 275  x time int = zeros(length(flux array{1,1})−1,1); xlen = length(x time int); for i = 1:(length(flux array{1,1})−1) x time int(i) = (x time(i+1)+x time(i))/2; end  149  D.2. Bead Capture Model Construction 276  % ... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−...  277 278  279  280  281  282  283  284 285 286 287 288 289  % Calculate linear distance from cell center to each bead ... center % ... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−... % because sortpos has been sorted the same way as the domain ... number is % sorted by comsol, the linearpos vector matches up with flux ... array (from % column 3 to column N+2. column 1 is the cell and column 2 is... all beads % combined in flux array and they aren't included here ... obviously linearpos = zeros(N,1); for i = 1:N linearpos(i) = sqrt(abs(sortpos(i,1)−cellpos(1))ˆ2+... abs(sortpos(i,2)−cellpos(2))ˆ2); end % ... −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−...  290 291 292  savnam = strcat('data out',num1,' ',num2,' ',num3,'.out'); neg space = −1;  293 294  295  296  save(savnam,'flux array sum','neg space','tot perc capt','... neg space',... 'tot intensity','neg space','mean intensity','neg space',... ... 'x time int','neg space','linearpos','−ASCII');  297 298  t1 = toc;  150  D.3. Bead Capture Data Analysis  D.3  1 2 3 4 5 6 7 8 9 10 11 12  13 14  Bead Capture Data Analysis  %% close all clear all clc %% flux sum = cell(3,4); tot perc capt = cell(3,4); tot intensity = cell(3,4); mean intensity = cell(3,4); x time int = cell(3,4); linear pos = cell(3,4); minfluxval = 5000; % minimum number of molecules to see light ... up r = 4.9/2; % bead radius in micron beadarea = 4*pi()*rˆ2;  15 16 17 18 19  20 21 22 23 24  25 26 27  28 29 30  31 32  data array temp = []; for i = 1:3 for j = 1:4 fname = strcat('data out',num2str(j),' ',num2str(i),' '... ,'1.out'); data array temp = dlmread(fname); [a,b]= size(data array temp); negpos = find(data array temp == −1); for k = 1:(negpos(1)−1) flux sum{i,j} = cat(1,flux sum{i,j},data array temp... (k,:)); end for k = (negpos(1)+1):(negpos(2)−1) tot perc capt{i,j} = cat(1,tot perc capt{i,j},... data array temp(k,1)); end for k = (negpos(2)+1):(negpos(3)−1) tot intensity{i,j} = cat(1,tot intensity{i,j},... data array temp(k,1)); end for k = (negpos(3)+1):(negpos(4)−1)  151  D.3. Bead Capture Data Analysis 33  34 35 36  37 38 39  40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57  58  59  60 61 62  mean intensity{i,j} = cat(1,mean intensity{i,j},... data array temp(k,1)); end for k = (negpos(4)+1):(negpos(5)−1) x time int{i,j} = cat(1,x time int{i,j},... data array temp(k,1)); end for k = (negpos(5)+1):a linear pos{i,j} = cat(1,linear pos{i,j},... data array temp(k,1)); end end end %% truemeanval = cell(3,4); truetotval = cell(3,4); horzFS len = cell(3,4); for i = 1:3 for j = 1:4 fluxS M = flux sum{i,j}; [s1,s2] = size(fluxS M); fluxS M = fluxS M(:,3:s2); [s1,s2] = size(fluxS M); for k = 1:s1 horzFS = fluxS M(k,:); fS1 = find(horzFS>=minfluxval); horzFS = horzFS(fS1); truetotval{i,j} = cat(2,truetotval{i,j},sum(horzFS)... ); truemeanval{i,j} = cat(2,truemeanval{i,j},sum(... horzFS)./(length(horzFS)*beadarea)); horzFS len{i,j} = cat(2,horzFS len{i,j},length(... horzFS)); end end end  63 64 65 66  figure for i = 1:4 subplot(2,2,i)  152  D.3. Bead Capture Data Analysis plot(x time int{1,1}./60,horzFS len{1,i},'r','linewidth',2) hold on plot(x time int{1,1}./60,horzFS len{2,i},'b','linewidth',2) plot(x time int{1,1}./60,horzFS len{3,i},'g','linewidth',2) ylim([0 140]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Included Beads','fontweight','bold') switch i case 1 titlen = 'Random Configuration'; case 2 titlen = 'Close Configuration'; case 3 titlen = 'Far Configuration'; case 4 titlen = 'Gradient Configuration'; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('60 Beads','100 Beads','140 Beads','location','... northwest'); grid on  67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89  90 91  end  92 93 94 95 96  97 98  99  100 101 102  figure for i = 1:4 subplot(2,2,i) plot(x time int{1,1}./60,truemeanval{1,i},'r','linewidth'... ,2) hold on plot(x time int{1,1}./60,truemeanval{2,i},'b','linewidth'... ,2) plot(x time int{1,1}./60,truemeanval{3,i},'g','linewidth'... ,2) ylim([0 1000]) xlim([0 10000/60]) set(gca,'fontweight','bold')  153  D.3. Bead Capture Data Analysis xlabel('Time (min)','fontweight','bold') ylabel('Measured Mean Intensity (molecules/umˆ2)','... fontweight','bold') switch i case 1 titlen = 'Random Configuration'; case 2 titlen = 'Close Configuration'; case 3 titlen = 'Far Configuration'; case 4 titlen = 'Gradient Configuration'; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('60 Beads','100 Beads','140 Beads','location','... northwest'); grid on  103 104  105 106 107 108 109 110 111 112 113 114 115 116 117 118  119 120  end  121 122 123 124 125 126 127 128 129 130 131 132 133  134 135 136 137 138 139  figure for i = 1:4 subplot(2,2,i) plot(x time int{1,1}./60,truetotval{1,i},'r','linewidth',2) hold on plot(x time int{1,1}./60,truetotval{2,i},'b','linewidth',2) plot(x time int{1,1}./60,truetotval{3,i},'g','linewidth',2) ylim([0 3.5e6]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Measured Total Intensity (molecules)','fontweight',... 'bold') switch i case 1 titlen = 'Random Configuration'; case 2 titlen = 'Close Configuration'; case 3  154  D.3. Bead Capture Data Analysis 140 141 142 143 144 145 146 147  148 149 150 151 152 153 154 155 156 157 158  159 160 161 162 163 164 165 166 167  168 169 170 171 172 173 174 175 176  titlen = 'Far Configuration'; case 4 titlen = 'Gradient Configuration'; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('60 Beads','100 Beads','140 Beads','location','... northwest'); grid on end figure for i = 1:4 subplot(2,2,i) valhold1 = []; valhold1(:,1) = truemeanval{1,i}; valhold1(:,2) = truemeanval{2,i}; valhold1(:,3) = truemeanval{3,i}; for j = 1:length(truemeanval{1,i}) percdiff1(j) = (max(valhold1(j,:))−min(valhold1(j,:)))... /((max(valhold1(j,:))+min(valhold1(j,:)))/2); end plot(x time int{1,1}./60,percdiff1.*100,'k','linewidth',2) hold on valhold2 = []; valhold2(:,1) = truetotval{1,i}; valhold2(:,2) = truetotval{2,i}; valhold2(:,3) = truetotval{3,i}; for j = 1:length(truetotval{1,i}) percdiff2(j) = (max(valhold2(j,:))−min(valhold2(j,:)))... /((max(valhold2(j,:))+min(valhold2(j,:)))/2); end plot(x time int{1,1}./60,percdiff2.*100,'b','linewidth',2) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Percent Difference','fontweight','bold') ylim([0,100]) xlim([0,10000/60]) switch i case 1  155  D.3. Bead Capture Data Analysis titlen case 2 titlen case 3 titlen case 4 titlen otherwise titlen  177 178 179 180 181 182 183 184  = 'Close Configuration'; = 'Far Configuration'; = 'Gradient Configuration';  = 'Error'; end title(titlen,'fontweight','bold') legend('Mean Intensity','Total Intensity') grid on  185 186 187 188 189 190  = 'Random Configuration';  end  191 192 193  %% floor total flux = dlmread('datfloor.txt');  194 195 196 197 198 199 200 201  202  203  204  for i = 1:(length(floor total flux)−1) floor total flux quad(i,1) = (floor total flux(i+1,1)−... floor total flux(i,1))*((floor total flux(i+1,2)+... floor total flux(i,2))/2); floor total flux quad(i,2) = (floor total flux(i+1,1)−... floor total flux(i,1))*325; floor total flux sum(i,1) = sum(floor total flux quad(:,1))... ; floor total flux sum(i,2) = sum(floor total flux quad(:,2))... ; floor time(i) = (floor total flux(i+1,1)+floor total flux(i... ,1))/2; end  205 206 207  vertlinex = ones(1000,1).*7200; vertliney = linspace(0,1,1000);  208 209  floor perc = floor total flux sum(:,1)./floor total flux sum... (:,2);  210 211 212  f = figure; set(f,'position',[130 100 1000 800])  156  D.3. Bead Capture Data Analysis 213 214 215 216 217  218  219  220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238  239 240 241 242 243 244 245 246  247  for i = 1:4 subplot(2,2,i) plot(floor time./60,floor perc,'−−k','linewidth',2) hold on plot(x time int{1,1}./60,tot perc capt{1,i},'r','linewidth'... ,2) plot(x time int{1,1}./60,tot perc capt{2,i},'b','linewidth'... ,2) plot(x time int{1,1}./60,tot perc capt{3,i},'g','linewidth'... ,2) ylim([0 1]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Percent Capture','fontweight','bold') switch i case 1 titlen = 'Random Configuration'; case 2 titlen = 'Close Configuration'; case 3 titlen = 'Far Configuration'; case 4 titlen = 'Gradient Configuration'; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('Floor','60 Beads','100 Beads','140 Beads','location... ','southeast'); grid on end f = figure; set(f,'position',[130 100 1000 800]) vertliney = linspace(0,800,1000); for i = 1:4 subplot(2,2,i) plot(x time int{1,1}./60,mean intensity{1,i},'r','linewidth... ',2) hold on  157  D.3. Bead Capture Data Analysis 248  249  250 251 252 253 254  255 256 257 258 259 260 261 262 263 264 265 266 267 268  269 270 271 272 273 274 275 276 277  278 279  280  plot(x time int{1,1}./60,mean intensity{2,i},'b','linewidth... ',2) plot(x time int{1,1}./60,mean intensity{3,i},'g','linewidth... ',2) ylim([0 800]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Mean Intensity (molecules/umˆ2)','fontweight','bold... ') switch i case 1 titlen = 'Mean Intensity − Random Configuration'; case 2 titlen = 'Mean Intensity − Close Configuration'; case 3 titlen = 'Mean Intensity − Far Configuration'; case 4 titlen = 'Mean Intensity − Gradient Configuration'; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('60 Beads','100 Beads','140 Beads','location','... northwest'); grid on end %% f = figure; set(f,'position',[130 100 1000 800]) vertliney = linspace(0,800,1000); for i = 1:4 subplot(2,2,i) plot(x time int{1,1}./60,mean intensity{1,i},'r','linewidth... ',2) hold on plot(x time int{1,1}./60,truemeanval{1,i},'−−r','linewidth'... ,2) plot(x time int{1,1}./60,mean intensity{2,i},'b','linewidth... ',2)  158  D.3. Bead Capture Data Analysis 281  282  283  284 285 286 287 288  289 290 291 292 293 294 295 296 297 298 299 300 301 302 303  304 305 306 307 308  309 310  311  plot(x time int{1,1}./60,truemeanval{2,i},'−−b','linewidth'... ,2) plot(x time int{1,1}./60,mean intensity{3,i},'g','linewidth... ',2) plot(x time int{1,1}./60,truemeanval{3,i},'−−g','... linewidth',2) ylim([0 1000]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Mean Intensity (molecules/umˆ2)','fontweight','bold... ') switch i case 1 titlen = 'Mean Intensity − Random Configuration'; case 2 titlen = 'Mean Intensity − Close Configuration'; case 3 titlen = 'Mean Intensity − Far Configuration'; case 4 titlen = 'Mean Intensity − Gradient Configuration'; otherwise titlen = 'Error'; end grid on title(titlen,'fontweight','bold'); legend('60 Beads','60 Beads − Measured','100 Beads','100 ... Beads − Measured','140 Beads','140 Beads − Measured','... location','northwest'); end %% for i = 1:3 figure plot(x time int{1,1}./60,mean intensity{i,1},'k','linewidth... ',2) hold on plot(x time int{1,1}./60,mean intensity{i,2},'b','linewidth... ',2) plot(x time int{1,1}./60,mean intensity{i,3},'r','linewidth... ',2)  159  D.3. Bead Capture Data Analysis plot(x time int{1,1}./60,mean intensity{i,4},'g','linewidth... ',2) grid on switch i case 1 titlen = 'Mean Intensity − 60 Beads'; case 2 titlen = 'Mean Intensity − 100 Beads'; case 3 titlen = 'Mean Intensity − 140 Beads'; otherwise titlen = 'Error'; end xlim([0 10000/60]) ylim([0 800]) xlabel('Time (min)','fontweight','bold') ylabel('Mean Intensity (molecules/umˆ2)','fontweight','bold... ') title(titlen,'fontweight','bold'); set(gca,'fontweight','bold') legend('Random','Close','Far','Gradient','location','... northwest')  312  313 314 315 316 317 318 319 320 321 322 323 324 325 326 327  328 329 330  331  end  332 333 334 335 336 337  338 339  340  341 342 343 344 345  f = figure; set(f,'position',[130 100 1000 800]) for i = 1:4 subplot(2,2,i) plot(x time int{1,1}./60,tot intensity{1,i},'r','linewidth'... ,2) hold on plot(x time int{1,1}./60,tot intensity{2,i},'b','linewidth'... ,2) plot(x time int{1,1}./60,tot intensity{3,i},'g','linewidth'... ,2) ylim([0 3.5e6]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold') ylabel('Total Intensity (molecules)','fontweight','bold')  160  D.3. Bead Capture Data Analysis 346 347 348 349 350 351 352 353 354  355 356 357 358 359  360 361 362 363 364 365 366 367  368 369  370  371  372  373  374 375 376 377  switch i case 1 titlen = 'Total Intensity − Random Configuration'; case 2 titlen = 'Total Intensity − Close Configuration'; case 3 titlen = 'Total Intensity − Far Configuration'; case 4 titlen = 'Total Intensity − Gradient Configuration'... ; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('60 Beads','100 Beads','140 Beads','location','... northwest'); grid on end %% f = figure; set(f,'position',[130 100 1000 800]) for i = 1:4 subplot(2,2,i) plot(x time int{1,1}./60,tot intensity{1,i},'r','linewidth'... ,2) hold on plot(x time int{1,1}./60,truetotval{1,i},'−−r','linewidth'... ,2) plot(x time int{1,1}./60,tot intensity{2,i},'b','linewidth'... ,2) plot(x time int{1,1}./60,truetotval{2,i},'−−b','linewidth'... ,2) plot(x time int{1,1}./60,tot intensity{3,i},'g','linewidth'... ,2) plot(x time int{1,1}./60,truetotval{3,i},'−−g','linewidth'... ,2) ylim([0 3.5e6]) xlim([0 10000/60]) set(gca,'fontweight','bold') xlabel('Time (min)','fontweight','bold')  161  D.3. Bead Capture Data Analysis ylabel('Total Intensity (molecules)','fontweight','bold') switch i case 1 titlen = 'Total Intensity − Random Configuration'; case 2 titlen = 'Total Intensity − Close Configuration'; case 3 titlen = 'Total Intensity − Far Configuration'; case 4 titlen = 'Total Intensity − Gradient Configuration'... ; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold'); legend('60 Beads','60 Beads − Measured','100 Beads','100 ... Beads − Measured','140 Beads','140 Beads − Measured','... location','northwest'); grid on  378 379 380 381 382 383 384 385 386 387  388 389 390 391 392  393 394  end  395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414  %% f = figure; set(f,'position',[130 100 1000 800]) horzlinex = linspace(0,180,1000); horzliney = ones(1000,1).*10000; for i = 1:4 subplot(2,2,i) fluxM = flux sum{1,i}; scatter(linear pos{1,i},fluxM(85,3:62),'r') hold on fluxM = flux sum{2,i}; scatter(linear pos{2,i},fluxM(85,3:102),'b') fluxM = flux sum{3,i}; scatter(linear pos{3,i},fluxM(85,3:142),'g') switch i case 1 titlen = 'Random Configuration'; case 2 titlen = 'Close Configuration';  162  D.3. Bead Capture Data Analysis case 3 titlen = 'Far Configuration'; case 4 titlen = 'Gradient Configuration'; otherwise titlen = 'Error';  415 416 417 418 419 420  end grid on title(titlen,'fontweight','bold') xlabel('Linear Distance from Cell (um)','fontweight','bold'... ) ylabel('Total Captured Molecules','fontweight','bold') xlim([0 200]) ylim([0 4.5e5]) legend('60 Beads','100 Beads','140 Beads') set(gca,'fontweight','bold')  421 422 423 424  425 426 427 428 429 430  end  431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453  %% for st = 1:3 switch st case 1 beadtit = '60'; case 2 beadtit = '100'; case 3 beadtit = '140'; otherwise end f = figure; set(f,'position',[130 100 1000 800]) for w = 1:4 fluxM = flux sum{st,w}; [s1,s2] = size(fluxM); for i = 3:s2 if fluxM(85,i)>=minfluxval over(i−2) = 1; else over(i−2) = 0; end  163  D.3. Bead Capture Data Analysis 454  end  455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481  482 483  484 485  486 487  488 489  linpos = linear pos{st,w}; linpos h = []; for i = 1:length(linpos) if over(i)==1; linpos h(i)=1; else linpos h(i)=0; end end linposz1 = find(linpos h>0); linposz2 = find(linpos h==0); linposp = linpos(linposz1); linposn = linpos(linposz2); subplot(2,2,w) x = [0:10:160]; hist(linposn,x); hold on hist(linposp,x); h = findobj(gca,'Type','patch'); set(h(1),'FaceColor','r'); set(h(2),'FaceColor','b'); xlim([0,200]) ylim([0,30]) switch w case 1 titlen = strcat({'Random Configuration − '},num2str... (beadtit),{' Beads'}); case 2 titlen = strcat({'Close Configuration − '},num2str(... beadtit),{' Beads'}); case 3 titlen = strcat({'Far Configuration − '},num2str(... beadtit),{' Beads'}); case 4 titlen = strcat({'Gradient Configuration − '},... num2str(beadtit),{' Beads'}); otherwise titlen = 'Error';  164  D.3. Bead Capture Data Analysis end title(titlen,'fontweight','bold') molval = strcat('Beads w/ <',num2str(minfluxval),'mol.'); molvalg = strcat('Beads w/ >',num2str(minfluxval),'mol.'); legend(molval,molvalg) xlabel('Linear Distance from Cell (um)','fontweight','bold'... ) ylabel('Number of Beads','fontweight','bold') grid on  490 491 492 493 494 495  496 497 498 499  end end  500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528  %% f = figure; set(f,'position',[130 100 1000 800]) maxb1 = 0; for w = 1:4 for u = 1:3 fluxM = flux sum{u,w}; [s1,s2] = size(fluxM); for i = 3:s2 if fluxM(85,i)>=minfluxval over(i−2) = 1; else over(i−2) = 0; end end linpos = linear pos{u,w}; linpos h = []; for i = 1:length(linpos) if over(i)==1; linpos h(i)=1; else linpos h(i)=0; end end numpos(u) = sum(linpos h); linposz1 = find(linpos h>0); linposp = linpos(linposz1); x = [0:10:160];  165  D.3. Bead Capture Data Analysis 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560  561  562 563 564 565 566  % % %  y1{1,u} = hist(linposp,x); yholdm = y1{1,u}; for st = 2:length(x) yholdm(st)=yholdm(st)/x(st); end maxyholdm = max(yholdm); scalef = 1/maxyholdm; yholdm = yholdm.*scalef; y1{1,u} = yholdm; end subplot(2,2,w) b1 = bar(x,[y1{1,1};y1{1,2};y1{1,3}]',1.5); maxb1hold = max(max(cell2mat(y1))); if maxb1hold>maxb1 maxb1 = maxb1hold; end b1h = get(b1,'Children'); set(b1h{1},'CDataMapping','direct'); set(b1h{2},'CDataMapping','direct'); set(b1h{3},'CDataMapping','direct'); mycolor=[0 1 0; 0 0 1; 1 0 0]; set(b1h{1}, 'CData',1); set(b1h{2}, 'CData',2); set(b1h{3}, 'CData',3); colormap(mycolor); num60 = strcat({'60 Beads − '},num2str(numpos(1))); num100 = strcat({'100 Beads − '},num2str(numpos(2))); num140 = strcat({'140 Beads − '},num2str(numpos(3))); legend(num60{1},num100{1},num140{1},'location','northeast') xlabel('Linear Distance from Cell (um)','fontweight','bold'... ) ylabb = strcat('Number of Beads with>',num2str(minfluxval),... 'mol. / rˆ2'); ylabel(ylabb,'fontweight','bold') switch w case 1 titlen = 'Random Configuration'; case 2  166  D.3. Bead Capture Data Analysis titlen = 'Close Configuration'; case 3 titlen = 'Far Configuration'; case 4 titlen = 'Gradient Configuration'; otherwise titlen = 'Error'; end title(titlen,'fontweight','bold') grid on  567 568 569 570 571 572 573 574 575 576 577  end  578 579 580 581 582 583 584 585 586 587  for i = 1:40 n5(i) = i*5; end i = 1; n5val=0; while maxb1>=n5val n5val = n5(i); i = i+1; end  588 589  maxb1 = n5val;  590 591 592 593 594 595 596 597 598 599 600 601 602  subplot(2,2,1) xlim([0 200]) ylim([0 .5]) subplot(2,2,2) xlim([0 200]) ylim([0 .5]) subplot(2,2,3) xlim([0 200]) ylim([0 .5]) subplot(2,2,4) xlim([0 200]) ylim([0 .5])  603 604 605 606  for k=1:17 figure(k); saveas(gcf,['fig',num2str(k)],'png');  167  D.3. Bead Capture Data Analysis 607  end  168  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0073985/manifest

Comment

Related Items