Medicine, Faculty of
Population and Public Health (SPPH), School of
DSpace
UBCV
McLaughlin, Angela
Master of Science - MSc
University of British Columbia
Background Despite the widespread availability of antiretroviral therapy, sub-populations are burdened by disproportionately high HIV transmission within both developed and un-developed countries. By identifying groups and areas where people are at a higher risk of HIV transmission, we can better target interventions to those communities. Phylogenetic analysis of HIV genetic sequences generated from resistance testing for people living with HIV can augment traditional epidemiological methods to quantify transmission. Each virus’ branching rate in the phylogenetic tree, referred to as its lineage-level diversification rate, is an estimate of its transmission rate. Methods The viral diversification rate was calculated for 36,271 partial pol HIV sequences representing 9,310 participants in the British Columbia Drug Treatment Program between May 1996 and March 2018. First, the relationships between patient-level diversification rate with age, viral load, health authority, risk behavior, ethnicity, and clinical background were explored using a generalized linear model. Subsequently, patient data was aggregated according to their census tract of residence in order to build a community-level zero-inflated negative binomial model to predict new HIV cases. Results At the patient level, high HIV diversification rate was positively associated with being younger, reporting injection drug use, having co-infection with hepatitis C virus, having a high most recent viral load, and residing within the Northern BC Health Authority, while having ever had AIDS and identifying as black were both significantly associated with lower diversification rates. At the census tract level, the sum of patients’ annual change in diversification rate improved the specificity, accuracy, and correlation of observed and predicted values of a model to predict new HIV cases relative to a model including only epidemiological and clinical patient data. Conclusions Analyzing the distribution of phylogenetic HIV diversification rate among people living with HIV in BC corroborated known risk factors for transmission and revealed areas at risk of new infections. This tool could be useful for public health surveillance and serves as proof of principle that phylogenetics could further inform the allocation of harm reduction services while still maintaining patient confidentiality.
APPLYING PHYLOGENETIC LINEAGE-LEVEL DIVERSIFICATION RATE TO INFER EPIDEMIOLOGICAL RISK FACTORS AND SPATIOTEMPORAL DISTRIBUTION OF HIV TRANSMISSION IN BRITISH COLUMBIA by Angela McLaughlin B.Sc. Microbiology and Immunology, The University of British Columbia, 2014 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE AND POSTDOCTORAL STUDIES (Population and Public Health) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2019 © Angela McLaughlin, 2019 ii The following individuals certify that they have read, and recommend to the Faculty of Graduate and Postdoctoral Studies for acceptance, a thesis/dissertation entitled: APPLYING PHYLOGENETIC LINEAGE-LEVEL DIVERSIFICATION RATE TO INFER EPIDEMIOLOGICAL RISK FACTORS AND SPATIOTEMPORAL DISTRIBUTION OF HIV TRANSMISSION IN BRITISH COLUMBIA submitted by Angela McLaughlin in partial fulfillment of the requirements for the degree of Master of Science in Population and Public Health Examining Committee: Jeffrey Joy Co-supervisor Jean Shoveller Co-supervisor Rodney Knight Supervisory Committee Member Jennifer Gardy Additional Examiner Additional Supervisory Committee Members: Richard Harrigan Supervisory Committee Member iii Abstract Background Despite the widespread availability of antiretroviral therapy, sub-populations are burdened by disproportionately high HIV transmission within both developed and un-developed countries. By identifying groups and areas where people are at a higher risk of HIV transmission, we can better target interventions to those communities. Phylogenetic analysis of HIV genetic sequences generated from resistance testing for people living with HIV can augment traditional epidemiological methods to quantify transmission. Each virus’ branching rate in the phylogenetic tree, referred to as its lineage-level diversification rate, is an estimate of its transmission rate. Methods The viral diversification rate was calculated for 36,271 partial pol HIV sequences representing 9,310 participants in the British Columbia Drug Treatment Program between May 1996 and March 2018. First, the relationships between patient-level diversification rate with age, viral load, health authority, risk behavior, ethnicity, and clinical background were explored using a generalized linear model. Subsequently, patient data was aggregated according to their census tract of residence in order to build a community-level zero-inflated negative binomial model to predict new HIV cases. Results At the patient level, high HIV diversification rate was positively associated with being younger, reporting injection drug use, having co-infection with hepatitis C virus, having a high most recent viral load, and residing within the Northern BC Health Authority, while having ever had AIDS and identifying as black were both significantly associated with lower diversification rates. At the census tract level, the sum of patients’ annual change in diversification rate improved the specificity, accuracy, and correlation of observed and predicted values of a model to predict new HIV cases relative to a model including only epidemiological and clinical patient data. Conclusions Analyzing the distribution of phylogenetic HIV diversification rate among people living with HIV in BC corroborated known risk factors for transmission and revealed areas at risk of new infections. This tool could be useful for public health surveillance and serves as proof of principle that phylogenetics could further inform the allocation of harm reduction services while still maintaining patient confidentiality. iv Lay Summary Despite widespread availability of antiretroviral therapy, there remain areas and demographic groups with disproportionately high HIV transmission. By identifying groups of people who are at an elevated risk of HIV transmission, we can target prevention and treatment services more efficiently. While estimating HIV transmission using contact tracing is hampered by uncertainty and biases, it is possible to estimate transmission activity by comparing patients’ viral genetic sequences to trace their viruses’ shared ancestry. When HIV infects a host, it branches into a distinct lineage in the evolutionary tree. This branching rate is the lineage-level diversification rate, which estimates viral transmission rate. The associations of sociodemographic and clinical factors with HIV diversification rate were explored. Then we built a model to predict new HIV infections by aggregating diversification rate and other clinical data by patient geography. These analyses could help to better allocate public health services to prevent and treat HIV. v Preface The idea of studying the HIV lineage-level diversification rate was prompted by Jeffrey Joy, who had previously applied this metric in other contexts. Analyzing diversification rate at an aggregate level to predict new HIV cases was my idea. The bioinformatics and statistical pipeline prior to and post obtaining diversification rates was written and executed by me. Jeff ran his code to calculate diversification rates. The data used was collected as part of the Drug Treatment Program at the BC Centre for Excellence in HIV/AIDS. This was the collective effort of the clinical lab team, data entry team, and software team at the Centre, as well as the physicians in contact with the patients. I wrote Chapter 2 and received edits from Jeffrey Joy. A version of this chapter may be submitted to support the publication of Chapter 3 and/or Chapter 4. I conducted the data analysis, interpretation, and writing for Chapter 3 and I received feedback from Jeffrey Joy and Richard Harrigan. Ethical approval for this study was granted by the University of British Columbia - Providence Health Care Research Ethics Board (H17-01812). This work has been presented at Conference for Retroviruses and Opportunistic Infections in Seattle in March 2019 and will also be presented at Canadian Conference on HIV/AIDS Research in Saskatoon SA in May 2019. A version of Chapter 3 will be submitted for publication in the New England Journal of Medicine. I conducted the data cleaning, analysis, interpretation, and writing for Chapter 4 and I received the most edits from Jeffrey Joy and some further suggestions for improvement from Richard Harrigan, Jeannie Shoveller, and Julio Montaner. The figure of the phylogenetic tree colored by lineage-level diversification rate was generated by Jeffrey Joy. An earlier version of the analysis was vi presented in poster and/or oral form at Conference for Retroviruses and Opportunistic Infections in Boston MA in March 2018, HIV Dynamics and Evolution in Leavenworth WA in April 2018, and at Canadian Conference on HIV/AIDS Research in Vancouver BC in April 2018. It was covered under my approved ethics, as granted by the University of British Columbia - Providence Health Care Research Ethics Board (H17-01812). A version of Chapter 4 will be submitted for publication in the New England Journal of Medicine. vii Table of Contents Abstract .......................................................................................................................................... iii Lay Summary ................................................................................................................................. iv Preface ..............................................................................................................................................v Table of Contents .......................................................................................................................... vii List of Tables ................................................................................................................................ xii List of Figures .............................................................................................................................. xiii List of Abbreviations .....................................................................................................................xv Acknowledgements ..................................................................................................................... xvii Chapter 1: Introduction ....................................................................................................................1 1.1 Epidemiology in the modern world ................................................................................ 1 1.1.1 The co-evolution of humans and pathogens ........................................................... 1 1.1.2 Phylogenetic approaches to prioritized health care ................................................ 2 1.2 The HIV epidemic at a global and local scale ................................................................ 3 1.2.1 The discovery of HIV as the causative agent of AIDS ........................................... 3 1.2.2 Origin and global distribution of HIV subtypes and groups ................................... 4 1.2.3 Global HIV epidemiology ....................................................................................... 5 1.2.4 The epidemiological history of HIV-1 in British Columbia ................................... 6 1.2.5 The BC Centre for Excellence in HIV/AIDS ......................................................... 9 1.3 The evolution of antiretroviral therapy ......................................................................... 10 1.3.1 Early antiretroviral therapies ................................................................................. 10 1.3.2 Highly active antiretroviral therapies .................................................................... 11 1.3.3 Treatment as prevention ........................................................................................ 12 viii 1.3.4 Pre-exposure prophylaxis ...................................................................................... 13 1.3.5 First- and second-line regimens in resource-limited settings ............................... 13 1.4 Challenges in treating and preventing HIV infections .................................................. 14 1.4.1 Imperfect treatment adherence .............................................................................. 14 1.4.2 Transmission of drug resistance and immune resistance ...................................... 15 1.4.3 Latent reservoirs .................................................................................................... 17 1.5 Phylogenetic analysis of HIV sequences ...................................................................... 18 1.5.1 Improvements in sequencing and computational technologies ............................ 18 1.5.2 An array of phylogenetics methods ...................................................................... 20 1.5.3 Transmission inference ......................................................................................... 22 1.5.4 Lineage-level diversification rate ......................................................................... 23 1.6 Ethics of HIV phylogenetics ......................................................................................... 24 1.6.1 Appropriate interpretations of phylogenetic data ................................................. 24 1.6.2 Criminalization of HIV non-disclosure ................................................................ 25 1.6.3 Researcher’s role in protecting patient confidentiality ......................................... 26 1.6.4 Researcher’s role in preventing ongoing transmission ......................................... 27 1.7 Targeted approaches to HIV eradication ...................................................................... 27 1.7.1 Attaining the UNAIDS 90-90-90 target for 2020 ................................................. 27 1.7.2 Research questions ................................................................................................ 28 1.7.3 Hypotheses ............................................................................................................ 29 1.7.4 Objectives ............................................................................................................. 29 1.7.5 Chapter guide ........................................................................................................ 29 Chapter 2: Exploring the roots of lineage-level diversification rate and the utility of phylogenetics in pathogen epidemiology .......................................................................................31 ix 2.1 Introduction ................................................................................................................... 31 2.2 History of the lineage-level diversification rate ............................................................ 32 2.2.1 Group level phylogenetic metrics ......................................................................... 33 2.2.2 Individual taxa phylogenetic metrics .................................................................... 33 2.2.3 Equal splits and the species-level lineage diversification rate .............................. 35 2.3 Applications of phylogenetics and lineage-level diversification rate in infectious disease epidemiology ................................................................................................................ 37 2.3.1 Comparing lineage-level diversification rate and genetic clustering .................... 38 2.3.2 Contemporary methods of estimating pathogen transmission rates ..................... 39 2.3.3 Potential applications of lineage-level diversification rates in epidemiology ...... 40 2.4 Conclusion .................................................................................................................... 42 Chapter 3: Epidemiological correlates of HIV lineage-level diversification rate ..........................43 3.1 Introduction ................................................................................................................... 43 3.2 Methods......................................................................................................................... 46 3.2.1 Study setting and participants ............................................................................... 46 3.2.2 Clinical, epidemiological, and demographic variables ......................................... 47 3.2.3 Sequences .............................................................................................................. 48 3.2.4 Phylogenetic inference .......................................................................................... 49 3.2.5 Diversification rate ................................................................................................ 49 3.2.6 Distribution of diversification rates across health authorities and age group ....... 50 3.2.7 Modelling .............................................................................................................. 50 3.3 Results ........................................................................................................................... 51 3.3.1 Patient population ................................................................................................. 51 3.3.2 Distribution of diversification rate by health authority and age category ............ 52 x 3.3.3 Epidemiological, clinical, and demographic corelates of diversification rate ...... 55 3.4 Discussion ..................................................................................................................... 57 3.5 Conclusion .................................................................................................................... 61 Chapter 4: Viral diversification rate identifies HIV transmission hot spots ..................................62 4.1 Introduction ................................................................................................................... 63 4.2 Methods......................................................................................................................... 64 4.2.1 Study setting and participants ............................................................................... 64 4.2.2 Phylogenetic inference .......................................................................................... 65 4.2.3 Diversification rate ................................................................................................ 66 4.2.4 Predictive modelling ............................................................................................. 66 4.3 Results ........................................................................................................................... 67 4.3.1 Study population ................................................................................................... 67 4.3.2 Diversification rate ................................................................................................ 69 4.3.3 Predictive modeling .............................................................................................. 70 4.4 Discussion ..................................................................................................................... 73 4.5 Conclusions ................................................................................................................... 76 Chapter 5: Conclusion ....................................................................................................................77 5.1 Overall analysis and integration of the research and conclusions ................................ 77 5.2 Re-visiting original thesis hypotheses .......................................................................... 78 5.3 Overall significance and contribution ........................................................................... 79 5.4 Strengths and limitations ............................................................................................... 80 5.5 Potential applications of the research findings ............................................................. 81 References ......................................................................................................................................83 Appendices .....................................................................................................................................94 xi Appendix A Supplementary material for Chapter 4: Viral diversification rate identifies HIV transmission hot spots ............................................................................................................... 94 A.1 Supplementary methods ............................................................................................ 94 Sequences .......................................................................................................... 94 Clinical, epidemiological, and demographic variables ..................................... 94 Diversification rate robustness .......................................................................... 95 Aggregated variables ........................................................................................ 96 Spatial analyses ................................................................................................. 97 Model training and testing ................................................................................ 98 A.2 Supplementary results and discussion ...................................................................... 99 A.2.1 Outcome validity ............................................................................................... 99 A.2.2 Diversification rate robustness and characterization ...................................... 100 A.2.3 Zero-inflated negative binomial model interpretation .................................... 102 A.2.4 Spatial analyses and interpretations ................................................................ 102 A.3 Supplementary data and shapefiles ......................................................................... 117 Statistics Canada census boundary files ......................................................... 117 Statistics Canada census tract data .................................................................. 117 Statistics Canada postal code conversion file ................................................. 118 Health authority boundary files ...................................................................... 118 xii List of Tables Table 1–1 Drug classes, approved drugs, and mechanisms. ......................................................... 11 Table 3–1 The study population included Drug Treatment Program participants in British Columbia with sequence data who were alive on March 29 2018. .............................................. 53 Table 4–1 Characteristics of the Drug Treatment Program cohort with sequences compared to the study population restricted to those who have reported their census tract of residence. .............. 68 Table 4–2 The final zero-inflated negative binomial (ZINB) predictive model. .......................... 71 Table A-1 A comparison of diversification rates calculated using different tree building strategies………………………………………………………………………………………..111 Table A-2 Pearson’s correlations between all evaluated potential aggregated predictors and the outcome. ...................................................................................................................................... 111 Table A-3 The goodness-of-fit of the full zero-inflated negative binomial (ZINB) model in comparison to competing models. .............................................................................................. 115 Table A-4 A comparison of the mean predictive power measures for the full zero-inflated negative binomial (ZINB) model to its reduced model without any diversification rate. .......... 114 Table A-5 The annual number of new HIV cases with sequences compared to those with census tract information available and estimated provincial incidence. ................................................ 116 xiii List of Figures Figure 3-1 Five regional health authorities in British Columbia. ................................................. 48 Figure 3-2 The association of age category and lineage-level diversification rate varies by health authority. ....................................................................................................................................... 54 Figure 3-3 A summary of the adjusted relative risks for patient attributes associated with HIV lineage-level diversification rates. ................................................................................................ 56 Figure 3-4 A summary of the adjusted relative risks for patient attributes associated with annual changes in HIV lineage-level diversification rate between 2018 and 2017. ................................. 57 Figure 4-1 A representative bootstrap approximate maximum likelihood phylogenetic tree with branches coloured by the lineage-level diversification rate. ......................................................... 69 Figure 4-2 Heat maps of the epicenter of the HIV epidemic in BC, downtown Vancouver, on odd study years between 2008 and 2017. ............................................................................................ 72 Figure 4-3 A comparison of the predictive power of the full ZINB and ZINB without phylogenetic measures (reduced ZINB). ...................................................................................... 73 Figure A-1 A comparison of annual estimates of new HIV cases between the BC Centre for Disease Control and this study……………………………………………………………... …..104 Figure A-2 The distribution of total annual new HIV cases within merged census tracts in BC across the study period. ............................................................................................................... 105 Figure A-3 The total number of tips in each time interval phylogenetic tree. ............................ 105 Figure A-4 The distribution of patients’ mean and % coefficient of variation (CV) for diversification rates calculated across 100 bootstrap approximate maximum likelihood trees. . 106 Figure A-5 Merged census tracts within British Columbia, Metro Vancouver, and downtown Vancouver. .................................................................................................................................. 107 xiv Figure A-6 Predictive power of the full zero-inflated negative binomial (ZINB) model compared to the reduced ZINB model over time for testing, training, and blind prediction data. .............. 108 Figure A-7 Workflow of data segregation, model training, testing, and blind predictions. ....... 109 xv List of Abbreviations AIDS = acquired immune deficiency syndrome aRR = adjusted relative risk ART = antiretroviral therapy BC = British Columbia BCCDC = British Columbia Centre for Disease Control BC-CfE = British Columbia Centre for Excellence in HIV/AIDS CDC = Centers for Disease Control CI = confidence interval CT = census tract DRM = drug resistance mutation DRV = drug resistant variants DTES = downtown eastside (of Vancouver) DTP = Drug Treatment Program ES = equal splits (metric) FN = false negative FP = false positive GIS = geographic information systems HAART = highly active antiretroviral therapy HCV = hepatitis C virus HIV = human immunodeficiency virus HLA = Human leukocyte-antigen IDU = injection drug user INSTI = integrase strand transfer inhibitor xvi KS = Kaposi’s sarcoma MMPP = Markov-modulated Poisson process MSM = men who have sex with men NGS = next generation sequencing NPV = negative predictive value NRTI = nucleoside reverse transcriptase inhibitor NNRTI = non-nucleoside reverse transcriptase inhibitor OR = odds ratio PCP = Pneumocystis carinii pneumonia PI = protease inhibitor PLHIV = people living with HIV PPV = positive predictive value PrEP = Pre-exposure prophylaxis RAxML = Randomized Axelerated Maximum Likelihood SIV = simian immunodeficiency virus TasP = Treatment as Prevention TN = true negative TP = true positive UNAIDS = Joint United Nations Programme on HIV/AIDS VL = plasma viral load (copies viral RNA/mL plasma) WHO = World Health Organization ZIP = zero-inflated Poisson ZINB = zero-inflated negative binomial xvii Acknowledgements Firstly, I would like to thank my supervisors, thesis committee members, and external examiner. Jeff - you have been a constant source of ideas, inspiration, and help. Thanks for sharing your time, expertise, and beer with me over the past two years. Jeannie - thank you for your mentorship, encouragement, and career advice. It means a lot to me to have your support. Richard - I can never thank you enough for taking me on at the Centre and putting me in touch with Jeff. Thank you for seeing my potential, for focusing my research vision, and for keeping me on my toes. Rod – thanks for your insights and kindness. Thanks also to my external examiner, Jenn, for helping to support my thesis completion amid a busy schedule and significant life changes. I am thankful to my funding providers for keeping me in the green. First and foremost, thank you to the BC Centre for Excellence in HIV/AIDS at St. Paul’s Hospital and Providence Health Care for employing me over the past two years. A huge thanks to the Canadian Institute for Health Research for selecting me for a Canada Graduate Scholarship in my first year of study. Thanks also to the UBC Faculty of Medicine for selecting me for the Dorothy Helmer award in 2018 for my research in the field of infectious disease. I also send my gratitude to the Health Officer’s Council of British Columbia for supporting me with the Dr. Ken Benson Memorial Bursary. To my colleagues in the Bioinformatics group and larger laboratory group at the BC Centre for Excellence in HIV/AIDS - thank you for working together to generate an unprecedented data set, not to mention helping people get personalized therapy based on your sequencing efforts. The xviii support I have received from my colleagues at the Centre has helped to teach me innumerable lessons about HIV sequencing, virology, programming, and evolution. To my professors at the UBC School of Population and Public Health – namely Dr. Amy Manges, Mike Marin, Dr. Eugenia Oviedo-Joekes, Dr. Trevor Dummer – thank you for teaching me how to think as an epidemiologist and prepare me for a career in public health. Additionally, I would like to thank my professors from the UBC Bioinformatics program – Dr. Paul Pavlidis and Dr. Sara Mostafavi - for expanding my mind to the power and limitations of bioinformatics. To my parents and siblings - thank you for helping me to believe in my own abilities, for fostering my curiosity with the natural world, and for supporting me through challenging times. Lastly, to my partner, Will - thank you for supporting me when I felt discouraged, over-whelmed, confused, and stressed (usually all at once). Thanks for listening to my gripes when things weren’t going as expected, as well as sharing in my excitement and enthusiasm when things went well. xix Dedication To the positive community, 1 Chapter 1: Introduction Clinical and research efforts to understand human immunodeficiency virus (HIV) and acquired immune deficiency syndrome (AIDS) have helped to catalyze immense changes in molecular biology, sequencing technology, epidemiological methods. The body of HIV and AIDS research produced in the past 36 years since the virus was discovered is concomitantly immense. Through this introduction, I draw upon only a fraction of the work done, focusing primarily upon the historical epidemiology of HIV globally and within British Columbia (BC), the challenges of treating or curing HIV infections, how genetic sequencing and phylogenetics have helped to illuminate transmission patterns, ethical considerations, and what work is needed to improve targeted HIV treatment and prevention services. Finally, I introduce my research questions, hypotheses, and objectives in the context of this body of research, as well as a chapter guide for the remainder of the thesis. 1.1 Epidemiology in the modern world 1.1.1 The co-evolution of humans and pathogens Our exponentially growing human population of over 7.6 billion has an expanding appetite for health care resources. While modern medicine breakthroughs such as antibiotics, vaccines, and hygienic practices have decreased the probability of dying from an infectious disease and extended the average life expectancy in many nations1, pathogens evolve alongside and within us. Although non-communicable diseases accounted for the large majority of deaths in high-income countries in 2016, communicable diseases - including lower respiratory infections, diarrheal diseases, HIV/AIDS, malaria, tuberculosis - still caused the majority of deaths in low-income countries1. The crude death rate of HIV/AIDS alone in low-income countries in 2016 was 2 over 40 per 100,000 people1. This is not to say that no progress has been made; the scale-up of antiretroviral therapy (ART) has driven a 48% decrease in HIV-related deaths globally, from a peak of 1.9 million in 2005 to 1 million in 20162. Continued efforts and funding will be required to further expand treatment and pre-exposure prophylaxis (PrEP) availability, as well as scale-up point of care testing and support other prevention and harm reduction services. The evolution and proliferation of drug resistance within many pathogens including HIV challenge us to develop new technologies and distribute resources effectively. The World Health Organization (WHO) identified HIV and antimicrobial resistance, respectively, as two of the top ten threats to global health in 20193. Sequencing clinical HIV isolates to screen for drug resistance mutations (DRMs) has become increasingly common in the past two decades and has resulted in unprecedentedly large datasets of viral sequences, often annotated with patient attributes. We now have the ability to track the evolution of HIV across space and time. One of the greatest challenges in modern epidemiology is to ensure that the distribution of resources among different socioeconomic and cultural groups is equitable and accessible. By monitoring the transmission patterns of infectious diseases, it is possible to estimate when, where, and how people were infected. 1.1.2 Phylogenetic approaches to prioritized health care With the advent of rapid genetic sequencing technologies and powerful computational methods, modeling pathogen transmission is increasingly a reality. Phylogenetics – the study of evolutionary relationships between biological entities – enables us to infer likely transmission networks based on the similarity of pathogens’ genetic sequences4,5. Although the phylogenetic methods applied to inform pathogen evolution inevitably carry some uncertainty, and opinions differ on which 3 models are most appropriate, they give us the ability to compare the evolutionary trajectory of biological entities, in this case viruses. Further, by identifying sub-populations or characteristics associated with heightened HIV transmission, phylogenetic HIV studies have resulted in real changes in the way prioritized treatment and prevention services are delivered6,7. 1.2 The HIV epidemic at a global and local scale 1.2.1 The discovery of HIV as the causative agent of AIDS In 1981, an unusually high incidence of Kaposi’s sarcoma (KS) and a rare opportunistic infection called Pneumocystis carinii pneumonia (PCP) among previously healthy homosexual men in Los Angeles followed by New York and San Francisco was brought to the attention of the Centers for Disease Control (CDC)8. The CDC assembled a task force, which determined that between June 1 1981 and March 24 1982 there were a total of 290 cases of KS, PCP, and other serious opportunistic infections with an overall 39.7% mortality rate9. Of these cases, 238 were homosexual or bisexual men, and of the remaining group, 56% were injection drug users9. Notably, 27 individuals were diagnosed with both KS and PCP, which had never before been reported9. That these rare diseases were primarily clustered among homosexual men of the same age group, racial distribution, and geographic areas strongly suggested that this was a single epidemic9. Immunosuppression was the common thread linking these two diseases. Differing hypotheses emerged about what the mysterious causative agent of the later coined acquired immune deficiency syndrome (AIDS) could be. Early case control studies and prospective cohort studies suggested that AIDS diagnosis in homosexual men was significantly associated with illicit drug use, history of mononucleosis, history of gonorrhea, frequency of anal receptive intercourse, age less than 30, and number of sexual partners in the preceding year10,11. Finally in 1983, a novel 4 retrovirus was isolated from the lymph node of a homosexual Caucasian man showing multiple lymphadenopathies and was determined to be the causative agent of AIDS12. Although several names were proposed for this novel virus including human T-lymphotropic virus type III (HTLV-III), lymphoadenopathy-associated virus (HAV), and AIDS-associated retrovirus (ARV), eventually an International Committee decided on its final nomenclature as human immunodeficiency virus (HIV)13. The more that we studied the genetic composition of HIV, the more we realized how diverse the virus really was. 1.2.2 Origin and global distribution of HIV subtypes and groups Evolutionary analyses have allowed researchers to trace back the ancestry of HIV to specific species of origin, as well as places and times of cross-over with reasonable certainty. We now know that at least eleven zoonotic events between non-human primates and primates gave rise to the HIV pandemic14. Among the primate lentiviruses family, HIV is most closely related to simian immunodeficiency viruses (SIV) found in non-human primates in Africa in which SIV infection usually does not result in a decreased CD4 count or manifest disease in their primate hosts15. HIV can be broadly characterized as either HIV-1, which is most closely related to SIV in chimpanzees16, or as HIV-2, which is most closely related to SIV in sooty mangabeys17. The HIV-1 lineages, Groups M, N, and O, are the result of at least 3 zoonotic events16. In contrast, HIV-2 strains are the result of at least 8 species jumping events, giving rise to Groups A-H, although only groups A and B are transmissible between humans and are primarily restricted to West Africa18. Though all HIV-1 strains are transmissible between humans, Group M is responsible for the majority of the HIV pandemic16. This lineage is estimated to have jumped from chimpanzees to humans around 1908 (95% confidence interval(CI), 1884-1924)19, likely in southeastern 5 Cameroon20. Within Group M, there are at least 9 subtypes – A, B, C, D, F, G, H, J and K - that exhibit diverse geographic distribution. For instance, whereas 98.3% of HIV infections in Southern Africa in 2004 were subtype C, only 0.45% of sampled HIV in North America were subtype C as 98.4% of infections were subtype B21. Subtype A dominates in Eastern Europe and central Asia (78.9%)21. Worldwide, subtype C accounted for 50% of infections in 2004, followed by 12% subtype A, 10% subtype B, 6% subtype G, 3% subtype D, and 13% circulating recombinant forms21. Looking at the distribution of HIV subtypes in BC in 2015, the large majority were subtype B (92.9%), with approximately 3.1% subtype C, 1.2% circulating recombinant form 01AE, 0.7% A1, 0.5% D, and 1.8% other subtypes22. 1.2.3 Global HIV epidemiology The Joint United Nations Programme on HIV/AIDS (UNAIDS) estimated the global HIV prevalence at the end of 2017 was 36.9 (95% CI, 31.1-43.9) million, 21.7 (95% CI, 19.9-22.6) million of whom were accessing treatment, and 1.8 (95% CI, 1.4-2.4) million of whom were new cases23. This represents a sizeable increase in treatment uptake from 2010, when only 8 (95% CI, 7.1-8.3) million PLHIV were on treatment23. Relative to the peak of the epidemic in 1996, new HIV diagnoses decreased by 47% by the end of 201723. In the context of the 90-90-90 targets set by UNAIDS for 2020, globally in 2017 75% of PLHIV were diagnosed, 79% of those diagnosed were accessing treatment, and 81% of those on treatment were virologically suppressed, resulting in a total of 47% of PLHIV being virally suppressed worldwide23. Inter-region epidemiological differences are vast. East and southern Africa experienced a 14% decline in new HIV infections between 2010 and 2015, which still accounted for 46% of the global total of new HIV infections24. Within sub-Saharan Africa in 2017, women aged 15-24 were 6 twice as likely to be living with HIV than men, and women who experienced sexual violence were 1.5 times more likely to become infected with HIV than other women23. Unfortunately, the eastern Europe and central Asia region experienced a 57% increase in new HIV infections between 2010 and 2015, and in 2015 only 21% (95% CI, 20-23%) of PLHIV in the region were on treatment24. Some of the largest global populations of injection drug users (IDUs) reside in eastern Europe and central Asia, with some countries estimated to have greater than 40% HIV prevalence among subpopulations of IDU25. In western and central Europe and North America, there were an estimated 2.2 (95% CI, 1.9-2.4) million people living with HIV in 2017, with 70,000 (95% CI, 57,000-84,000) new diagnoses23. Globally, the risk of acquiring HIV in 2017 was 27 times higher for men who have sex with men (MSM), 23 times higher for IDU, 13 times higher for female sex workers, and 13 times higher for transgender women23. Additionally, prison populations were reported to have 5 times the HIV prevalence than the general population in North America26. Every region and population affected by HIV tells a different epidemiological story, with varying dominant transmission routes, genetic variants, health care service availability, and community support. In light of this diversity, different geographies demand specific strategies that target the right groups for improved connection to treatment and prevention services. 1.2.4 The epidemiological history of HIV-1 in British Columbia The HIV epidemic in British Columbia (BC) has been dynamic and complex. In the early 1980’s, the majority of new cases were among MSM, particularly among those with high numbers of partners and who engaged in frequent receptive anal intercourse27. By the 1990’s, the epidemic had shifted into the IDU community primarily located in the downtown eastside (DTES) of 7 Vancouver28. In conjunction with the viral prevalence among IDU, female sex workers were also deemed to be at a significantly higher risk of HIV seroconversion29. With the provision of free highly active antiretroviral therapy (HAART) in concert with drug resistance testing and other harm reduction infrastructure, HIV incidence has steadily decreased in the province30,31. BC has been very progressive in spearheading harm reduction programs to decrease the incidence of HIV. As early as 1988, a needle exchange program was initiated in the DTES of Vancouver29. Later in 2013, the first medically supervised injection facility in North America, Insite, was opened in Vancouver. Insite was conservatively estimated to prevent 5-6 HIV infections per year32, in addition to serving as a first point of contact with HIV+ individuals who may otherwise be disconnected with care33. The notion of Treatment as Prevention (TasP) was empirically validated in BC30. Much of the overall decrease in HIV incidence in BC is attributable to a decrease in the number of new diagnoses among IDU; by 2016, only 6.6% of new HIV diagnoses were among IDU31. Although great strides have been made to decrease overall HIV incidence, the rate of new diagnoses for MSM in BC has not declined as rapidly as for other groups31. In 2016, MSM accounted for over 60% of new diagnoses in BC31. Young MSM are the greatest concern – the 1990-1999 and 1980-1989 birth cohorts both saw an increase in the number of HIV diagnoses between 2015 and 201631. This trend is not unique to BC. One potential reason may be that young individuals are more likely to engage in risky sexual behavior with partners who they met through applications like GRINDR, a social networking application popular among young MSM for sexual partnering34. Additionally, these age cohorts have grown up in an era where HIV has become a 8 manageable chronic disease, reducing fear of acquiring HIV and potentially leading to poorer precautionary measures. First Nations people also continue to be disproportionately affected by the HIV epidemic. The rate of new diagnoses among First Nations people in 2016 was 15.2 cases per 100,000, which exceeded the overall rate of new diagnoses of 5.1 per 100,00031. Although the percent of all new HIV cases in BC who identified as First Nations has decreased to 9.1% in 2016, 28.9% of all new female diagnoses in BC are First Nations women31. Furthermore, the rate of reported AIDS cases among First Nations people has consistently exceeded the overall rate in the BC population between 2007 and 2016, and the rate of AIDS case reports among both male and female First Nations people have increased between 2015 and 201631. There are many reasons that First Nations and other Indigenous people experience inequities in health access, which makes this a complex, multifaceted problem which presents special challenges. In particular, colonization, intergenerational trauma, residential schools, loss of language and culture, geographic isolation, discrimination, and poverty have all been cited as potential barriers to accessing HIV services for Indigenous people31. It is imperative that health care providers foster trust with Indigenous people and strive to deliver care that is consistent with their needs and world views. Although models project that BC will meet the UNAIDS 90-90-90 targets for 202035, particular groups including young MSM and female First Nations people remain at a high risk of contracting HIV. We should strive to provide more tailored services to these groups in close proximity to where they reside and make efforts to engage them with treatment and prevention. 9 1.2.5 The BC Centre for Excellence in HIV/AIDS Since opening in 1992, the British Columbia Centre for Excellence in HIV/AIDS (BC-CfE) at St. Paul’s Hospital has become an internationally recognized, world leading, multidisciplinary institution that is highly engaged in both clinical and research activities. The laboratory at the BC-CfE operates an extensive next-generation sequencing platform for drug resistance testing, providing this service to citizens from across BC, as well as other provinces and nations. In addition to being a world leader in providing personalized antiretroviral therapy, the BC-CfE was instrumental in showing the effectiveness TasP30 and supporting the UNAIDS 90-90-90 targets for 202035. The objectives of the BC-CfE36 are to “to improve the health of British Columbians through the development, ongoing monitoring and dissemination of comprehensive investigative and treatment programs for viral diseases”; “to develop new methods of patient care through vigorous clinical trials on prevention, investigation, treatment and management of HIV infection and HIV related diseases”; “to develop rational and cost-effective investigational and therapeutic protocols and programs for the prevention and treatment of HIV and related viral diseases”; “to provide education and other support programs to health care professionals”; “to participate in on-going surveillance and epidemiology programs; and “to maintain effective cost-benefit analysis programs related to the investigation and treatment of viral diseases”. The Drug Treatment Program (DTP) - the data source for the analyses at hand - is one of several cohorts based out of the BC-CfE. It maintains an extensive patient data set comprised of HIV genetic sequences, plasma viral load measurements, treatment history, epidemiological, and sociodemographic data. 10 1.3 The evolution of antiretroviral therapy 1.3.1 Early antiretroviral therapies Amid a mounting disease burden in the 1980’s, the first antiretroviral therapy (ART) was approved for clinical use. Zidovudine, also known as azidothymidine (AZT), is a thymidine analogue within the broader class of nucleoside reverse transcriptase inhibitors (NRTIs). AZT was shown to block viral replication37, prolong survival, decrease the frequency and severity of opportunistic infections, and transiently improve CD4 counts38. Unfortunately, the initially recommended dose of 1,500 mg was toxic for many, resulting in adverse reactions such as severe anemia, headaches, weight loss, and nausea, all of which were lessened at lower dosages38. In addition to a number of other NRTI therapies approved thereafter, entirely new drug classes were also developed, including non-nucleoside reverse transcriptase inhibitors (NNRTIs), protease inhibitors (PIs), fusion inhibitors, and integrase strand transfer inhibitors (INSTIs). Table 1–1 summarizes currently available drugs and brief mechanisms of each class39. 11 Table 1–1 Drug classes, approved drugs, and mechanisms. Drug class Approved drugs Mechanism Nucleos(t)ide analogue reverse transcriptase inhibitor (NRTI) Zidovudine Abacavir Didanosine Emtricitabine Lamuvudine Stavudine Tenofovir A nucleoside or nucleotide analogue lacking a 3’-hydroxl group is incorporated in the new strand of viral DNA by HIV reverse transcriptase, blocking further elongation. Non-nucleoside reverse transcriptase inhibitor (NNRTI) Delaviridine Efavirenz Nevirapine Etravirine Rilpivirine Doravirine A molecule that directly binds to hydrophobic NNRTI binding pocket on reverse transcriptase, preventing it from reverse transcribing viral RNA into DNA. Protease inhibitor (PI) Amprenavir Atazanavir Darunavir Fosamprenavir Indinavir Lopinavir Nelfinavir Ritonavir Saquinavir Tipranavir PIs bind HIV protease, blocking proteolytic cleavage of gag and gag-pol polyproteins necessary for HIV to become a mature and infectious virus. Fusion (entry) inhibitor Enfuvirtide Maraviroc Interfere with binding, fusion, and/or entry of HIV into host cells. Enfurvitide binds to gp41, while maraviroc prevents the interaction of HIV with the CCR5 co-receptor. Integrase strand transfer inhibitor (INSTI) Raltegravir Elvitegravir Dolutegravir Bictegravir INSTIs block HIV integrase, preventing HIV from integrating its viral DNA into the host genome and preventing viral replication. 1.3.2 Highly active antiretroviral therapies As early as 1996, the benefits of combined therapy were documented. For instance, a randomized clinical trial of 2,467 HIV-1-infected individuals showed that treatment with AZT and an additional NRTI slowed the progression of AIDS and improved survival relative to treatment with AZT alone40. An international consensus of therapeutic guidelines was established as a result, which recommended to “monitor plasma viral load and CD4 count; initiate treatment before immunodeficiency becomes apparent; aim to reduce plasma viral load as much as possible for as long as possible; use combinations of at least two drugs (which maximize antiviral effect, minimize 12 cross-resistance, and avoid overlapping toxicity); and change to a new regimen if plasma viral load rebounds despite continued therapy” 41. Further, it was found three drug regimens resulted in lower plasma viral loads than two drug regimens amongst antiretroviral-naive individuals42. Combined antiretroviral therapy, also commonly referred to as highly active antiretroviral therapy (HAART) has become the standard of care in developed nations where there are adequate resources for drug resistance testing to inform personalized combinations of ARTs. Modern HAART regimens usually include two NRTIs and one NRTI or PI39; alternatively, a combination of INSTIs or entry inhibitors with RT inhibitors (NRTI or NNRTI) and PIs can be used39. HAART has led to sustained decreases in morbidity, mortality, and transmission among people living with HIV (PLHIV)30. Administering HAART to PLHIV suppresses viral replication (plasma HIV-1 RNA levels <50 copies/mL) leading to decreases in the concentration of HIV-1 RNA in plasma (viral load)42 and other biological fluids including semen43, vaginal fluids44, and anal mucosa45. Individuals undergoing treatment consequently have a decreased probability of transmitting HIV vertically46, sexually47, and by sharing needles48. 1.3.3 Treatment as prevention HAART is effective not only to reduce the morbidity and mortality associated with HIV, but also to decrease transmission at the population level30,48. While this seems obvious to us now, it was not immediately evident whether wide availability of treatment would lead to an increase in risky behaviors that would negate any potential benefits49. Adherence to HAART has been demonstrated to decrease the concentration of HIV-1 RNA in biological fluids, resulting in a decreased risk of vertical transmission, sexual transmission, and blood-borne transmission48. There are no documented cases of individuals with undetectable viral load transmitting HIV sexually 13 (undetectable=untransmissible)50. Montaner et al. were among the first to identify that an increase in HAART coverage was associated with decreases in community viral load and the number of new HIV diagnoses per year51. From an economic perspective, HAART provides an additional cost savings when considering its protective effect on HIV transmission. For example, the provision of HAART to 3,693 PLHIV in BC in 2005 cost $49 million and was estimated to have averted 400 new infections in BC, which represented an averted lifetime expenditure of $96.4 million49. 1.3.4 Pre-exposure prophylaxis Starting in 2012, the WHO recommended oral pre-exposure prophylaxis (PrEP) containing tenofovir disoproxil fumarate for HIV serodiscordant couples and MSM52. In a systematic review of eighteen studies evaluating evidence for the use of oral PrEP containing tenofovir disoproxil fumarate in populations at elevated risk for HIV, Fonner et al. found that PrEP was significantly protective against HIV infection across populations, resulted in few safety risks, and did not result in changes in sexual risk behaviour53. Although there were more cases of drug-resistant HIV infection among PrEP users who initiated PrEP while acutely HIV-infected, the incidence of acquiring drug-resistant HIV during PrEP use was low53. 1.3.5 First- and second-line regimens in resource-limited settings In resource-limited countries lacking the funding, personnel, and/or infrastructure to conduct drug resistance testing, WHO recommends a simplified first-line and second-line regimen for adults54. The most commonly prescribed regimens are tenofovir (an NRTI) with lamivudine (an NRTI) or emtricitabine (an NRTI) with efavirenz (an NNRTI), and if those are contraindicated, zidovudine 14 (an NRTI) with lamivudine (an NRTI) and nevirapine (an NNRTI) is recommended54. If a patient fails their first-line regimen as measured by virologic failure, a second-line regimen is then recommended based on the patients’ first-line regimen, with the expectation that the previously unused NRTI will be effective. However, accumulated drug resistance mutations (DRMs) can compromise the effectiveness of the second-line NRTI backbone, particularly if virological failure is not immediately detected. If a patient’s first-line regimen contained zidovudine, the odds of having a compromised second-line regimen are 9.9 times higher than if they had a first line regimen containing tenofovir54. There are concerns that the ramp up of first-line regimens in developed nations will lead to parallel increases in the prevalence of DRMs, particularly for NRTIs39. Monitoring viral load and transmission of drug resistant variants within populations in the developing world is critical to contain the epidemic, yet unfortunately does not occur at an adequately expansive scale. 1.4 Challenges in treating and preventing HIV infections 1.4.1 Imperfect treatment adherence In 2014, it was estimated that 94% of PLHIV in BC who were on treatment were virally suppressed35. While this is reason to be encouraged, there remains a population on treatment within whom HIV actively replicates. Imperfect adherence to HAART is a large contributor to the virally unsuppressed population. Adherence to HAART on a daily basis can be challenging for many PLHIV, particularly those who experience adverse side effects, unstable housing, or other traumatic life disturbances. They may also experience distrust of the medical system, trauma with healthcare providers, or large travel distances to access care. 15 There is a positive correlation between HAART adherence and sustained virologic suppression (plasma viral load <50 copies/mL), as well as CD4 cell count increase of at least 145 cells/µL55. Additionally, those with the lowest adherence had an increased odds of mortality and emergence of drug resistance55,56. Further, non-adherent patients on efavirenz-based NNRTI therapies were 3.61 times more likely to die than adherent patients and those on boosted PI-based regimens were 3.25 times more likely to die than adherent patients56. Imperfect adherence to HAART is a strong driver of virological rebound, emergence of drug resistance, AIDS, and death56. When a patient interrupts their treatment, latent viruses inaccessible to ARTs come out of dormancy and cause viral rebound. As the viruses in the reservoir are an archive of the within-host evolution that may have been under selective drug pressure57,58, it follows that ART-experienced patients experiencing viral rebound should have a higher likelihood of harboring drug resistant variants than ART-inexperienced patients. PLHIV with relatively high but imperfect prescription refill percentages (80-90%) were found to be at an elevated risk of harbouring drug resistant HIV variants than those who had low prescription refill percentage (0-20%)59. 1.4.2 Transmission of drug resistance and immune resistance Drug resistant variants (DRVs) are selected at the within-host level and if they maintain replicative fitness, are transmitted between hosts. While patients are screened for DRVs when they are connected with care or after virologic failure in many countries, this is not always the case in the developing world where the majority of the global HIV burden resides. International efforts to support resistance testing for countries in need have been made by labs like the BC-CfE, but 16 realistically there are not enough resources and funding to do this for the 37 million PLHIV in the world. In a study of transmitted drug resistance conducted by Rhee et al. across 111 countries, between 2000 and 2013 the odds of NNRTI-associated transmitted drug resistance increased in North America, Europe, Latin America, and sub-Saharan Africa60. In sub-Saharan Africa, this time frame coincided with national ART scale-up, suggesting that first-line regimens were promoting the transmission of NNRTI resistance60. In low and middle income countries, the long term success of first-line regimens could be compromised by the low genetic barrier of resistance for NNRTIs60. Four DRMs accounted for >80% of all NNRTI transmitted resistance in all regions and subtypes, and sixteen DRMs accounted for >69% of all NRTI transmitted drug resistance60. In BC in 2005, those who initiated NNRTI-containing HAART had a 1.84 times higher risk of developing drug resistance59. Additionally, 25% of antiretroviral-naïve adults initiating HAART in BC harboured some level of transmitted drug resistance and the proportion of participants with DRMs increased during the administration of HAART59. In addition to evolving in response to drug pressure, HIV variants evolve in response to selective pressure from the immune system. Within a multicenter cohort of >1500 chronically infected treatment naïve PLHIV, 282 codons in gag, pol, and nef commonly mutated in response to Human Leukocyte-Antigen (HLA)-associated immune pressure61. Fully characterizing HLA-associated polymorphisms can help to inform common immune escape variants and vaccine discovery. Recently, HLA-adapted HIV variants were found to be significantly enriched in 17 transmission clusters in Saskatchewan, a province in which approximately 80% of HIV diagnoses are from people of Indigenous decent62. 1.4.3 Latent reservoirs The persistence of drug resistant and immune resistant HIV variants is dependent on the establishment of the latent reservoir, which themselves pose the greatest barrier to an effective HIV cure. Even in the presence of ART, HIV persists in long-lived reservoirs of latently infected cells, evading detection by the immune system and ART activity. Although CD4+ T cells are the primary reservoir, other cell types including macrophages and hematopoetic progenitor cells have been shown to harbour latent HIV – it is unclear how much they contribute to viral rebound63. The latent reservoirs represent an archive of an individuals’ within-host viral diversity over time58. This supports the feasibility of the ‘store and retrieve’ hypothesis that viruses stored in the latent reservoir early in infection preferentially transmit and establish new infections64. Although by definition, latency is the persistence of integrated HIV that is transcriptionally silent but still able to replicate, recent improvements in assays to distinguish intact from defective proviruses have revealed that the majority of HIV DNA integrated in host genomes in latent reservoirs is not fully intact or able to replicate57. Eliminating replication competent virus from reservoirs could be possible with a shock and kill approach in which a latency reversing agent shocks HIV out of latency and then antibodies, T-cell vaccines, or immunotherapies are used to promote the turnover of reactivated cells63. An effective cure or vaccine depends on our understanding of the composition and dynamics of latent reservoirs. 18 1.5 Phylogenetic analysis of HIV sequences Studying the evolutionary relatedness of HIV genetic sequences has fundamentally contributed to our understanding of the dynamics of between-host transmission and within-host evolution. Here, I briefly review current sequencing and computational technologies, common phylogenetic inference methods, how trees can estimate transmission networks, and the lineage-level diversification rate, which is reviewed in greater depth in Chapter 2. 1.5.1 Improvements in sequencing and computational technologies Sanger sequencing, also known as the chain termination method, was originally developed in 1977 and is one of the earliest forms of genetic sequencing. The template DNA is mixed with labelled deoxynucleotides and dideoxynucleotides and a DNA polymerase. As elongation occurs, some strands incorporate a dideoxynucleotide, which terminates elongation. In contemporary Sanger methodologies, strands are then separated by length on an electrophoresis gel and the terminating base for each strand is identified by laser excitation and spectral emission analysis. Sanger sequencing has accuracy up to 99.99%, but is limited in its throughput and cost65. The development of next generation sequencing (NGS) technologies, such as the Illumina platform, increased the throughput and decreased the cost per sample of sequencing65. While there are many NGS platforms that be either categorized as sequencing by ligation or sequencing by synthesis, they all share the characteristic of being massively parallelized to collect information from millions of reaction centres simultaneously. Accuracy is a little bit lower than with Sanger sequencing, but is still acceptably high at greater than 99%65. Read lengths are considerably shorter than Sanger sequencing, ranging from 75 to 700bp65. Although target amplification leads to high 19 sequencing depth, de novo assembly of short read lengths can result in non-contiguous consensus sequences, particularly in the presence of repetitive sequences and high % GC content65. Third generation sequencing technologies, such as the Oxford Nanopore MinION or the Pacific Biosystems single molecule real-time platform, represent the leading edge of sequencing. The nanopore technology fundamentally differs from its predecessors in that the detection of nucleotides is based on changes in electrical charge when each molecule passes through a nanopore embedded in a membrane. This technology supports ultra-long read lengths, in certain documented cases greater than 2 Mb66. Longer reads permit simpler assembly of a contiguous consensus sequence, particularly for sequences with a high % GC content or long repeats. However, sequencing depth is markedly lower and the accuracy varies depending on the bioinformatics pipeline used67. Pipelines that combine Illumina and nanopore approaches by first using the long nanopore reads for assembly and then polishing with short Illumina reads can reap the benefits of both read depth and length66. The nanopore platform is well-suited to outbreak response in the field as it is highly portable, light (100 g), does not require target amplification if there is an adequate quantity of starting genomic DNA67, and is supported by products for automated library preparation (Oxford Nanopore VolTRAX) and data analysis68. The volume of genomic data being generated as a result of advances in sequencing technologies is unprecedented and demands novel solutions for data storage, bioinformatics pipelines, analysis, and visualization. 20 1.5.2 An array of phylogenetics methods In a phylogenetic tree, tips or leaves represent biological entities - such as species or viral lineages – that can be traced back along shared branches all the way down to the root, which is the most recent common ancestor shared by all the tips. Where branches diverge, so do the evolutionary trajectories of the tips that emerge from those branches. Branching points are referred to as internal nodes and represent the inferred ancestors common to those branches. Reconstructing a phylogeny is akin to comparing differences between sampled leaves in order to infer the structure of the underlying tree that gave rise to those leaves. Most often, the sampled leaves only represent a fraction of the total leaves in the tree, thus an inferred tree is inherently an approximation of the true tree that gave rise to the sampled leaves. From that reconstructed tree, we can make further inferences about the rate of branching for different leaves, clusters of similar leaves, and even the time and location of when the tree’s root gave rise to the tree. While many types of character data such as morphology or functionality can be used to compare biological entities (leaves), the comparison of their genetic sequences is particularly powerful to reconstruct a molecular phylogeny to estimate their shared ancestry. A well-informed phylogeny permits us to look into the past to estimate the timing and genetic composition of lineages’ shared ancestry. While the implications of this are vast in evolutionary biology and pathogen epidemiology, it is not as simple as plug and play. Assumptions must be made, and models chosen. Ultimately, tree building is a balance of computational efficiency and thorough searching of the tree space. The collection of all possible trees, known as the tree space69, for most genomic alignments is too vast to search exhaustively. Tree-building algorithms apply different strategies to heuristically sample this space as efficiently as possible. 21 Frequentist (maximum likelihood) and Bayesian (conditional) approaches vary in their fundamental philosophies and results. Bayesian methods rely on prior assumptions about the underlying model giving rise to the data70. The BEAST software provides a Bayesian statistical framework for parameter estimation and hypothesis testing of many complementary evolutionary models70. Although specifying priors can in some cases be a benefit and improve model calibration, if no prior is obvious in the literature, then the researcher is burdened with ensuring that prior does not bias the sampled posterior distribution. Model misspecifications of tree priors in a molecular clock analysis, for example, can upwardly bias the clock rate71. BEAST uses a Markov chain Monte Carlo algorithm to sample from the posterior distribution for which the researcher must specify appropriate operators that propose new model parameters, otherwise it will not efficiently sample the distribution72. Current implementations of Bayesian inference are too computationally expensive for large alignments of tens of thousands of sequences. Maximum likelihood methods differ from Bayesian methods because they allow the data alone to inform the best model, as priors are inextricably linked to biases and may limit the generalizability of the results. Randomized Axelerated Maximum Likelihood (RAxML) is a popular program for building phylogenies for relatively large datasets using maximum likelihood73. Its primary strength is that it returns trees with good likelihood scores, although it is limited for large alignments74. The time for tree inference increases exponentially with the number of sequences and bootstrapping slows it down even further (100 bootstraps will take 100 times longer)74. For prohibitively long or wide alignments, approximate maximum likelihood implemented in FastTree 2 is a good alternative to RAxML74. It applies a minimum evolution subtree pruning and regrafting strategy to find a topology that minimizes the amount of evolution (or cumulative branch lengths). Although FastTree is less accurate than RAxML, most of the 22 disagreements are at nodes with poor support74. What FastTree lacks in accuracy, it makes up for with computational speed; it is at least 100 times faster than RAxML for alignments larger than 500 sequences74. 1.5.3 Transmission inference A phylogeny containing sequences from a well-sampled epidemic provides an estimate of the transmission network5. In one of the few instances where the time and direction of transmission within a cluster was known, the reconstructed trees were good estimates of the true phylogeny4. Identifying clusters of closely related HIV sequences has become a common method to detect subpopulations exposed to higher transmission rates than the general population6,75. There are a multitude of methods applied to build clusters from phylogenetic trees and results can vary widely depending on the method8,9. For instance, the chosen genetic distance threshold for cluster detection has a strong effect on which clusters are recovered from the true transmission network76,77. Some methods link two individuals in a network if their sequences are less than 1.5% distant (Tamura-Nei 93 distance, TN93)7,78, while others use a patristic distance (tree path length from one tip to another) cut-off based on the 95th percentile of within-host pairwise distances6. Another group of methods first identify subtrees with minimum bootstrap support values to then evaluate pairwise distances within each subtree79,80. Cluster detection is pragmatic for public health efforts as individuals within large or quickly growing clusters can be connected to prioritized care including partner services, provision of PrEP to partners, and linkage to treatment6,7. Additionally, it is possible to evaluate which risk groups are more likely to be in a cluster as a proxy for their probability of being involved in ongoing transmission78. However, cluster detection has been criticized for detecting differences in sampling, rather than differences in transmission in 23 subpopulations81. Furthermore, evaluating whether a patient is within a cluster or not is an over-simplification of the transmission dynamics within clusters. Not all individuals within a cluster are contributing towards new cases to the same extent; there is a spectrum of transmission activity within clusters. Other methods attempt to quantify the evolution of pathogen transmission rates across phylogenies using either a Markov-modulated Poisson process82 or a multi-type birth–death branching model process83, however they rely on wide assumptions about the evolution of rates along the tree. Some methods go so far as trying to connect transmission pairs and infer directionality84. Beyond simulated data, there has been limited testing of how robust and sensitive these methods are using a comprehensively sampled epidemic in humans where timing and direction of transmission are known81. In a widely-spanning epidemic, there are always unsampled hosts in a population and un-sampled variants within a host. These limitations should always be considered when interpreting phylogeny in the context of transmission. 1.5.4 Lineage-level diversification rate In this research, we validate and apply lineage-level diversification rate to estimate between-host transmission rates. The origins, calculation, and utility of lineage-level diversification rate in infectious disease epidemiology are described in greater detail in Chapter 2. Regardless of the method applied, any use of patients’ HIV sequence data to make inferences about possible transmission must consider the ethics of protecting patients’ privacy and confidentiality while also maximizing public benefit. 24 1.6 Ethics of HIV phylogenetics As next-generation sequencing becomes faster and more affordable, the potential to inform public health intervention using pathogens’ genetic sequences is increasingly feasible. With great power comes great responsibility. Here, I describe the appropriate interpretation of phylogenetics in HIV surveillance, review the landscape of HIV non-disclosure criminalization in Canada, and tease apart the balanced roles of the researcher to protect the interests of people living with HIV and those at risk of infection. A comprehensive review of the ethical considerations of HIV phylogenetic research was recently published by Coltart et al.85. 1.6.1 Appropriate interpretations of phylogenetic data While phylogenetics can be used to identify sub-populations exposed to higher transmission rates than the general population, individual transmission events cannot be inferred with certainty in all cases4. It is challenging to decipher transmission directionality and to dismiss the possibility that a third party infected both individuals or served as an intermediary86. In instances where multiple sequences are available from a putative transmission pair and the donor’s sequences are paraphyletic (the most recent common ancestor of all the sequences from an individual is shared with monophyletic subgroups from another individual) while the recipient’s sequences are polyphyletic (derived from more than one common ancestor), then directionality can be inferred with high certainty and intermediary links can be excluded87. There will always be some uncertainty in a phylogeny and thus phylogenetic reconstruction of transmission events in the context of forensic analyses and criminal law should be interpreted as such. Further, if the analysis was not conducted properly, a phylogeny could be misinterpreted 25 to make two sequences seem more similar than they are by including inappropriate controls88. Any prosecution case should draw upon a combination of different clinical and epidemiological evidence, including individuals’ sexual health histories, in addition to phylogenetic inference88. 1.6.2 Criminalization of HIV non-disclosure Concern has previously been raised that phylogenetic data could be misused in criminal cases of HIV non-disclosure89. Under current Canadian criminal law, PLHIV can be charged and prosecuted if they do not disclose their HIV positive status to their partner(s) before engaging in sexual activity that poses a realistic possibility of transmission, regardless of viral load90. The Supreme Court of Canada does not realize the oxymoron therein, as there have been no instances of PLHIV with an undetectable viral load transmitting the virus sexually50. Of further concern is that criminal charges for HIV non-disclosure can be laid even if HIV is not transmitted90. Until the inherent flaws in this legal guidance are rectified, we will continue to have to treat HIV differently than other notifiable diseases and efforts must be made to obscure the granularity of HIV phylogenetic analyses. There is legal precedent of phylogenetic information being used as evidence in the court of law to convict PLHV of not only HIV non-disclosure, but also sexual assault and in the most extreme cases, homicide88. One of the first times phylogenetics was applied within forensics was the investigation of HIV transmission from a dentist with AIDS in Florida in 1990 to dental surgery patients, following the HIV diagnosis of a patient with no other transmission risk exposure91. The dentist died during the investigation and no criminal charges were ever laid, however it is widely accepted that the dentist infected at least six patients with HIV. More recently, two index cases in the US were identified on the basis of paraphyly of sources sequences with respect to recipients92. 26 The criminalization of HIV non-disclosure has negative unintended consequences. Firstly, fear of prosecution disincentivizes diagnostic testing for HIV infection. If an individual does not know that they are positive, they are not liable90. It also promotes stigma around being HIV positive, which could lead to less disclosure overall to partners, counsellors, and frontline workers. Further, it counteracts the ‘undetectable=untransmissible’ message that is important to teach both the positive community and greater population. Most importantly, there is no evidence that criminalization of HIV non-disclosure has helped to prevent new infections90 1.6.3 Researcher’s role in protecting patient confidentiality Notwithstanding potential criminal charges of HIV non-disclosure, any leak of private patient information into the public domain could lead to social stigma and embarrassment for PLHIV. Fear of inappropriate data distribution could lead some individuals to refuse consent to have their data used for research purposes. Any breach of private information could deteriorate the trust of the HIV+ community with their health care providers and greater research community. In regard to publishing sensitive patient data, identifying patient data should never be included. To further safeguard against any type of privacy breach, patient data should always be anonymized in researchers’ datasets. Data access must be strictly regulated to personnel who are trained and knowledgeable about the data in a greater context, as well as to facilities and work stations with appropriate controls to limit any privacy breaches. Another way to ensure that individuals’ data remains anonymous to the public eye is by analyzing clusters, not individuals. Be they genetic clusters or geographic clusters, actionable results can still be generated without pointing a finger at any specific individual as being a potential transmitter. 27 1.6.4 Researcher’s role in preventing ongoing transmission As public health researchers, it is our ethical obligation both to protect the privacy of PLHIV as well as to prevent transmission among the population at risk. This includes assessing and mitigating potential public harm, and equally assessing the potential for public benefit. For instance, if phylogenetic reconstruction of a local epidemic reveals actively growing genetic clusters or the potential of a super-transmitter, it is our responsibility to communicate these findings to relevant public health authorities in a discrete way that protects that patient from any legal liability while also preventing further transmission. By the same token, an ethical researcher should be an ally for the positive community by educating others about HIV as a manageable chronic condition, the principle of ‘undetectable = untransmissible’, and effective harm reduction measures. 1.7 Targeted approaches to HIV eradication 1.7.1 Attaining the UNAIDS 90-90-90 target for 2020 With the aided influence of the BC Centre for Excellence in HIV/AIDS, the Joint United Nations Programme on HIV/AIDS (UNAIDS) set global targets towards the eradication of HIV. By 2020, nations should aim to have 90% of PLHIV diagnosed, 90% of those who are diagnosed on treatment, and 90% of those who are on treatment virologically suppressed. By 2030, nations should aim to have 95-95-9524,35. In BC, it was estimated in 2014 that 83% of PLHIV were diagnosed, 78% of those who are diagnosed were on treatment, and 84% of those on treatment were virologically suppressed35. 28 Connecting the sub-groups in the population who are not diagnosed, on treatment, or virologically suppressed is critical to achieve the UNAIDS goals and control HIV/AIDS once and for all. The heterogeneity of the pandemic requires strategies targeted specifically to each region, considering its dominant transmission routes and genetic variants. Phylogenetic analyses of HIV genetic sequences permit us to trace the evolution of HIV and connect hosts by the similarity of their viruses. When used in conjunction with clinical and epidemiological patient data, phylogenetics can be used to infer transmission events and timing with some certainty. Through my research, I sought to illuminate risk factors and predict spatiotemporal transmission patterns by analyzing the HIV lineage-level diversification rate of PLHIV in BC alongside their epidemiological and clinical data. By identifying transmission risk factors and hotspots, it is my hope that additional treatment and prevention services could be recommended within communities most at risk. 1.7.2 Research questions • What is the history and utility of the lineage-level diversification rate? • Is the HIV lineage-level diversification rate significantly correlated with known (and unknown) epidemiological, clinical, and demographic risk factors for HIV transmission? • When aggregated by patients’ geography of residence, is the HIV lineage-level diversification predictive of where and how many new cases arise in each community? • Is the predictive model of new HIV cases consistent with observed new data? 29 1.7.3 Hypotheses • HIV lineage-level diversification rate of PLHIV in BC is significantly associated with known transmission risk factors such as injection drug use and being young, in addition to previously unknown putative risk factors. • Geographically-aggregated HIV lineage-level diversification rate is a significant predictor of where new HIV cases arise in the future. • The predictive model trained and tested using the 2008-2017 data is predictive of ongoing data patterns through 2018. 1.7.4 Objectives • Review the history and potential usefulness of the lineage-level diversification rate. • Conduct a cross-sectional observational study within the Drug Treatment Program cohort in order to assess the association of epidemiological, clinical, and demographic cohort characteristics with HIV transmission rate, as estimated by viral lineage-level diversification rate. • Investigate whether patients’ geographically-aggregated viral diversification rate is predictive of communities where new HIV cases arise across British Columbia. • Validate the predictive model on new data for 2018. 1.7.5 Chapter guide In Chapter 2, I reviewed the history of lineage-level diversification rate in evolutionary biology and conservation biology and explored its potential applications in infectious disease epidemiology. In Chapter 3, I quantified the association of patients’ clinical, epidemiological, and 30 demographic attributes with lineage-level diversification rate. In Chapter 4, patients’ diversification rate, epidemiological, clinical, and demographics data were aggregated by their area of residence to predict new HIV cases across census tracts in British Columbia. Finally, in Chapter 5 I synthesized my results and compared my findings to existing research. 31 Chapter 2: Exploring the roots of lineage-level diversification rate and the utility of phylogenetics in pathogen epidemiology Angela McLaughlin1,2 and Jeffrey B. Joy1,3 1BC Centre for Excellence in HIV/AIDS, St. Paul’s Hospital, Vancouver, BC 2School of Population and Public Health, University of British Columbia Vancouver, BC 3Department of Medicine, University of British Columbia, Vancouver, BC Molecular phylogenies have permitted improved methods of quantifying the extent and rate of genetic divergence between biological entities. Comparing rates of genetic diversification between entities is useful in the context of evolutionary ecology, conservation biology, and pathogen epidemiology. Lineage-level diversification rate is a tip-level metric that reflects the splitting frequency of a lineage across its inferred ancestry, weighted more strongly towards the tip. Here, we review the origin of this metric and its utility in infectious disease epidemiology in comparison with other phylogenetic methods of inferring transmission. 2.1 Introduction Complementing epidemiology with phylogenetics contributes to our understanding of pathogen evolution and can lead to novel ways of delivering health services. In the case of the HIV epidemic, phylogenetic analyses of viral genetic sequences isolated from patients and annotated with their clinical and epidemiological information can shed light on active transmission clusters22,6,78. Individuals who are associated with a transmission cluster can thus be identified for enhanced connection to care, partner services, and treatment adherence follow-up6. Although cluster detection is becoming a common public health tool to identify high transmission groups to prioritize interventions for individuals within clusters, there is no agreed upon method of clustering. Results differ by method4,5 and common cluster detection methods have been criticized for primarily detecting differences in sampling rather than differences in transmission81 . 32 Studying the viral lineage-level diversification rate across a population overcomes this sensitivity to cluster definition. Rather, each tip in a viral phylogenetic tree is assigned a lineage-level lineage diversification rate, which is an estimate of the splitting frequency of that virus along the entire tree, weighted by proximity to the present94. Since transmission of HIV to a new host is equivalent to the formation of a new lineage in an evolutionary tree, diversification rates inferred from viral phylogenetic trees are a proxy for between-host transmission rates95. An important characteristic of lineage-level diversification rate is that the contribution made by each branch is weighted by its proximity to the tip94. Consequently, recent splitting events that are characterized by short terminal branch lengths increase the diversification rate more than events that occurred longer ago, which leads to a more current estimate of active transmission. Interestingly, lineage-level diversification rates have not been extensively applied in infectious disease epidemiology literature95. Diversification rates and other measures of phylogenetic distinctness have been mostly applied in conservation biology. Here, we explore the history of phylogenetic estimates of species distinctness and splitting frequency in order to better understand the origin, meaning, and alternatives to the lineage-level lineage diversification rate. Finally, the utility of the lineage-level diversification rate in epidemiology and other transmission estimation methods are discussed in the context of HIV and other infectious diseases. 2.2 History of the lineage-level diversification rate Inferring evolutionary relationships between taxa has long been utilized to inform the prioritized conservation of evolutionarily distinct species or habitats. There are two approaches that have been taken: the first quantifies the phylogenetic diversity held by a group of taxa, whereas the second focuses on elucidating the evolutionary distinctiveness of individual taxa. 33 2.2.1 Group level phylogenetic metrics Looking first at the group level phylogenetic metrics, phylogenetic diversity was one of the first to be incorporated into a conservation framework by Daniel Faith in 199296. This clade level metric is calculated by summing the appropriate branch lengths from the cladogram for a given set of taxa. While this metric is useful to compare relative diversity between sets of taxa, it is a limited to comparing sets of taxa because the phylogenetic diversity of each species on a cladogram, phylogram, or ultrametric phylogenetic tree is the same since all tips are the same distance from the root97. Magallon and Sanderson extended the concept of cladistic diversification rate to assess whether angiosperms’ relatively high species richness compared to other plants was attributable to different rates of speciation and extinction98. Clade-level diversification rate in this model was defined as the difference between the speciation rate and the extinction rate. Since there is often a lack of information about a clade’s diversity over time, the cladistic diversification rate was estimated using clades’ current age and species diversity. However, by definition, this clade level metric lacked the ability to assign rates to individual tree tips, regardless of the comparison being made. 2.2.2 Individual taxa phylogenetic metrics An alternative school of thought was that phylogenetic distinctness could be quantified at an individual species (or other taxa) level. Early conceptualization of species-level evolutionary distinctiveness, such as the Vane-Wright metric, only took tree topology into account, without considering branch lengths99. The Vane-Wright metric was partially informative because species that are more evolutionarily distinct have fewer nodes between the tip and root, but it suffers from poor resolution since it doesn’t incorporate branch lengths100. However, a subsequent extension, the May metric, both incorporated branch lengths and allowed consideration of polytomies101. An 34 alternative early estimate of species’ evolutionary isolation is the pendent edge, which is equivalent to the terminal branch length leading to the species tip102. While pendent edge values inherently ignore the internal tree topology, they have the advantage of being an absolute measure that can be compared between any taxa. Notably, pendent edge is sensitive to alternative species designations97. Contemporary estimations of species-level evolutionary distinctiveness consider both topology and internal branch lengths; they include the fair proportion metric, a quadratic entropy-based metric, and equal splits103,100,104. These node-based measures of evolutionary distinctiveness were previously reviewed by Pavoine et al.100, and they were further compared by Redding et al. in terms of how well they captured phylogenetic diversity97, so they will only be briefly compared here. The fair proportion metric, also known as the evolutionary distinctiveness score103, divides all internal branch lengths equally between all species below it regardless of tree structure103. Fair proportion has been further extended to incorporate species abundance, as per the abundance-weighted evolutionary distinctiveness105. One advantage of the fair proportion index is that it is independent of clade size in trees with greater than 100 species, making scores comparable between unrelated taxa103. In contrast, quadratic entropy is an explicitly pairwise measure between taxa that was first devised in 1982 by Rao106, and has since been expanded by Pavoine et al. as a quadratic entropy-based index100. Quadratic entropy is a measure of the dissimilarity between two species, calculated as the sum of branch lengths connecting them to a common ancestor. In a study comparing how well each index captured phylogenetic diversity, quadratic entropy, fair 35 proportion, Vane-Wright, and pendent edge metrics were consistently outperformed by equal splits97. 2.2.3 Equal splits and the species-level lineage diversification rate Redding and Mooers pioneered the equal splits (ES) metric, which is a measure of the evolutionary isolation of individual tips on a phylogenetic tree that takes branch length and topology into account104. The phylogenetic distance (in units of time or substitutions per site) of each edge is divided equally among its daughter edges until the entire tree has been sub-divided (Equation 1). The ES for a given tip is the sum of sub-divided branch lengths from root to tip, weighted by each branch’s proximity to the tip, such that the sum of all ES measures across the tree equals the sum of all the branch lengths107,97. When distinct species are represented in the tree tips, as for Redding and Mooers, the ES approximates the genetic distinctness that each species carries. Thus, a species that does not share much of its genetic history with other species in the tree will have a greater ES score than other species. The pendent edge metric described above is incorporated into the ES measure and for Yule trees it can contribute over half the ES score97. Notably, Redding and Mooers made the case for using time-scaled phylogenies to estimate ES because it is comparable across taxa and is intuitive for policy makers and laymen. 𝐸𝑆 =$ 𝑙&∏ (𝑑(𝑘 − 1) − 1)&./01&/0 𝐸𝑆2 = $ 𝑙&2&4056&/0 Equation 1. Equal splits (ES), as first devised by Redding and Mooers104, where j is the internal node on direct path from i to root (r), lj is the edge length from internal node j to j -1, and d(k) is the degree (three for bifurcation) at node k. Equation 2. Equal splits (ES) for a rooted bifurcating tree, as described by Jetz et al.94, where Ni is the number of edges n the path from species i and the root, and lj is the length of the edge j, with j=1 being the pendant edge leading to the species and j=Ni being the edge nearest the root. 36 1/𝐸𝑆2 = 𝐷𝑅2 = ($ 𝑙&2&4056&/0 )40 Jetz et al. extended the usage of the ES metric to investigate the species-level variation in diversification rate across a tree of all extant avian species94. In addition to quantifying ES across all birds, which for a bifurcating rooted tree is shown in Equation 2, they also analyzed the distribution of its reciprocal, the lineage-level lineage diversification rate, as in Equation 394. This lineage-level metric was shown to be complementary to the clade level net diversification rate described by Magallon and Sanderson. Diversification rate can be interpreted as the integrated splitting rate of a lineage from root to tip for any biological entity represented on the tree. An important feature of lineage-level diversification rate is that it is weighted more strongly towards the tip, emphasizing recent evolutionary events. Jetz et al. reported that diversification rate followed an approximately log-normal distribution94. They found that certain species contributed more towards rapid clade diversification and these rapid radiations were well distributed across the tree. Further, lineage-level diversification rates were spatially heterogeneous, such that Eastern hemispheres displayed lower average diversification rates than Western hemispheres. However, they found no significant difference in average diversification rates across latitudes. Areas characterized by lots of climatic fluctuations tended to have higher diversification rates and larger relative prevalence of top diversification rate species. This is consistent with another finding that non-pelagic bird species had significantly greater diversification rates than birds that generally inhabited the mainland. Their results supported that the main drivers of diversification are multiple and ongoing, individually constrained adaptive radiations. Equation 3. The species-level lineage diversification rate (DR) for each tip i on a rooted, bifurcating tree is the reciprocal sum of its sub-divided branch lengths lj across each branch j from the tip (j=1) to the root (j=Ni). 37 Conservation biologists have validated that lineage-level diversification rate and other phylogenetic metrics of evolutionary distinctness are informative and robust. By additionally analyzing the spatial distribution of phylogenetic distinctness, these metrics can further inform prioritization of conservation resources. Phylogenetic measures of evolutionary divergence have also been conceptualized and studied on large geographic scales. For instance, spatial mapping of imperiled phylogenetic endemism (IPE) among lemurs in Madagascar was used to identify areas containing spatially rare phylogenetic diversity additionally weighted by the probability of extinction108. Various phylogenetic and conservation parameters were mapped onto equi-rectangularly projected raster cells and then rank correlation coefficients and count sum rank tests were used to assess correlations and auto-correlation, respectively. In other studies, georeferenced branch extents have been applied to map phylogenetic endemism109,110. 2.3 Applications of phylogenetics and lineage-level diversification rate in infectious disease epidemiology The lineage-level diversification rate has only had limited application in infectious disease epidemiology, restricted so far to HIV95. When looking at a tree where each tip represents a virus sampled from an infected population, the phylogeny approximates the transmission network, and thus lineage-level diversification rates approximate transmission rates between hosts. Importantly, by weighing branching events closer to the tip (present day), transmission events are emphasized while still also considering the evolutionary history of the lineage. 38 Here, we review the predominant phylogenetic methods currently applied in HIV epidemiology, compare them conceptually to lineage-level diversification rate, and explore potential applications of this under-utilized metric in pathogen epidemiology. 2.3.1 Comparing lineage-level diversification rate and genetic clustering Genetic clustering is applied in infectious disease epidemiology to identify groups of closely related sequences, which should in theory correspond to subgroups with higher transmission activity than the general population. In public health practice, clustering of HIV sequences in particular has proven to be a useful tool to identify groups of individuals with augmented transmission levels22,78,7,111. These individuals can then be prioritized for enhanced partner contact services and pre-exposure prophylaxis dispensation, prioritized treatment follow-up, and other harm reduction services. There is little consensus in the field about what cluster detection method should be applied. Methods vary widely in regard to the genetic distance they use, how they establish a distance threshold, whether they incorporate bootstrap support values in their identification of clusters, the minimum number of individuals in a cluster, and whether they use subtrees to calculate clusters. These methods have been reviewed extensively in other literature6. The clusters identified by different methods are variable6,112 and there is a distinct lack of comprehensively sampled epidemics to empirically validate these methods. Comparing the methods’ performance in simulated epidemics, however, has indicated that all cluster detection methods are more sensitive to detecting differences in sampling than in transmission6. 39 While genetic clustering of rapidly-evolving pathogen sequences has certainly contributed towards our understanding of transmission dynamics, they have limitations. Lineage-level diversification rate used in conjunction with clusters could further leverage and validate existing prioritized care for individuals within large and rapidly growing clusters. 2.3.2 Contemporary methods of estimating pathogen transmission rates Alternative methods have been developed to infer disease transmission rates from phylogenetic trees, which differ from the lineage-level diversification rate in their assumptions about the evolution of rates along the tree. McCloskey and Poon, for instance, used a Markov-modulated Poisson process (MMPP) to model the evolution of transmission rates along the phylogeny to infer transmission clusters under the premise that changes in transmission rate are indicative of localized outbreaks82. While the MMPP model had higher specificity and sensitivity than competing clustering methods on simulated datasets, their results were framed around a series of assumptions. Namely, the MMPP model assumes that the rates at which lineages transition between branching rates are constant over time and that branching events themselves occur at a constant rate over time82. In light of the cumulative effects of epidemiogical, within-host, and sampling dynamics on the phylogeny, these assumptions seem hard to justify. An alternative method was previously introduced by Stadler and Bonhoeffer whereby transmission dynamics in structured populations were estimated using a multi-type birth-death branching model process, where each host is assigned to a sub-population83. The results were highly dependent on the assumption that a given hosts’ transmission rate, death rate, and recovery rate are common to their subpopulation. Their methods have been further expanded upon to infer clusters based on changes in transmission rates 77. Lineage-level diversification rate fundamentally 40 differs from these methodologies as there are no assumptions about the evolution of rates and the estimation of transmission rates is an end in itself instead of a way to identify clusters. 2.3.3 Potential applications of lineage-level diversification rates in epidemiology At the level of individual patients, lineage-level diversification rate can be used to investigate patient characteristics associated with augmented transmission rates in generalized linear models95. Similar studies have been conducted using clusters to approximate a higher transmission rate78,111, however by using cluster membership as their binary outcome, these analyses lacked resolution about individual differences in transmission rate within clusters and were themselves subject to the same biases as cluster detection itself. Our analyses have so far indicated that individuals who are treatment-experienced had a significantly lower risk of having a high diversification rate, supporting the notion of Treatment as Prevention95. Caution must be exercised when interpreting the lineage-level diversification rate as a proxy for transmission rate. An individual could have a high viral lineage-level diversification rate without themselves being an active transmitter if they are in close phylogenetic proximity to one or more individuals who are active transmitters. Through their shared viral ancestry, they share historic lineage diversification, and although the metric is weighted more strongly towards branching events close to the tip, the possibility exists for a lineage to have a high diversification rate driven by early rapid and extensive branching events, rather than recent ones. Researcher discretion is also required when reporting lineage-level diversification rates to ensure patients’ privacy and confidentiality is maintained. As not disclosing HIV status during any potential transmission event (even if an individual is virally suppressed and thus unable to 41 transmit) is criminalized in many nations, any patient data that could be identifying must be protected and anonymized wherever possible. Having the HIV+ community’s trust that their private data will be discretely utilized and for their benefit is critical for ongoing research in this field. Until legislation changes regarding HIV transmission, it is the responsibility of the researcher to protect patients’ sensitive data and clearly interpret lineage-level diversification rates as estimates of cumulative viral transmission activity. Aggregating patient data by their geographic area of residence is one way of protecting patient confidentiality, while maintaining resolution regarding differences in transmission activity. Spatial mapping of HIV lineage-level diversification rates can provide a way to focus prioritized care in areas where active transmission is occurring between hosts. In Chapter 4, we demonstrated that geographically-aggregated HIV lineage-level diversification rate, in conjunction with epidemiological and clinical parameters, was predictive of where and how many new HIV cases would arise in small geographic areas in the future. There also exists potential for lineage-level diversification rate to be used to explore within-host viral evolution. Examining the rate at which viral lineages diverge across different viral reservoirs over time could be a novel way of investigating which reservoirs produce the greatest variety of viral variants and the underlying characteristics of the reservoirs that predict diversification rates. Application to other rapidly-evolving pathogens besides HIV has yet to be explored, but certainly presents a novel opportunity, as well. We could learn more about the transmission dynamics of Zika, Ebola, and hepatitis C virus by investigating either the spatiotemporal 42 distribution of lineage-level diversification rates or looking at patient characteristics related to transmission activity. 2.4 Conclusion Following its origins in conservation biology and evolutionary ecology, lineage-level diversification rate illuminates closely related viruses or organisms, considering the entire inferred path of evolution. In doing so, it provides a robust estimation of between-host transmission rates when the tree is built using clinically-sampled pathogen sequences with broad implications for prioritized connection with care, health care resource distribution, and our understanding of pathogen evolution. 43 Chapter 3: Epidemiological correlates of HIV lineage-level diversification rate Angela McLaughlin1,2 and Jeffrey B. Joy1,3 1BC Centre for Excellence in HIV/AIDS, St. Paul’s Hospital, Vancouver, BC 2School of Population and Public Health, University of British Columbia Vancouver, BC 3Division of Infectious Diseases, Department of Medicine, University of British Columbia, Vancouver, BC Identifying HIV transmission risk factors can inform the prioritization of health care services. Although phylogenetic clustering of HIV sequences is routinely performed to evaluate characteristics associated with transmission, clusters membership over-simplifies the range of transmission activity across a population. We introduce an alternative phylogenetic method to investigate HIV transmission risk factors in British Columbia (BC), Canada, based on viral lineage-level diversification rates. For 8,103 people living with HIV (PLHIV) in BC in March 2018, we recovered the baseline HIV protease and RT sequences from the BC Drug Treatment Program database. We inferred 100 bootstrap approximate maximum likelihood phylogenetic trees and for each tip, we calculated its mean diversification rate as an estimate of its transmission rate. Patient attributes significantly associated with high viral transmission were evaluated using a gamma generalized linear model of cross-sectional 2018 data. Having a high HIV diversification rate was positively associated with being male, younger, using injection drugs, ever having hepatitis C virus, having a high recent viral load, and residing within the Northern BC Health authority (p<0.001). In contrast, having ever had AIDS and identifying as black were both significantly associated with lower diversification rates (p<0.001). By identifying risk factors associated with HIV transmission using the viral diversification rate among PLHIV, we can confidently recommend prioritized provision of treatment and prevention services for key groups. Additionally, these analyses highlight that HIV diversification rates can illuminate differences in individuals’ transmission activity, regardless of phylogenetic cluster membership. 3.1 Introduction Identifying risk factors for HIV transmission is critical for determining how to most effectively allocate limited public health resources for prevention and treatment. In British Columbia (BC) in 2014, it was estimated that 83% of people living with HIV (PLHIV) were diagnosed, 78% of those 44 diagnosed were on treatment, and 84% of those on treatment were virally suppressed35. While projections suggested that we are on track to meet the Joint United Nations Programme on HIV/AIDS (UNAIDS) 90-90-90 targets for 2020 and 95-95-95 targets for 203035, particular groups in the population remain disproportionately affected by HIV/AIDS. For instance, sub-populations in BC including women, young people, and injection drug users (IDU) were less likely to take up antiretroviral therapy (ART) and to achieve viral suppression35. It is additionally of interest to evaluate which subgroups are at higher risk for HIV transmission in order to better connect them to prioritized care. Quantifying HIV transmission dynamics by reconstructing social or sexual networks via traditional epidemiological methods such as partner tracing are limited by the virus’ long incubation period and low transmission rate per contact5. Conducting phylogenetic analyses of viral genetic sequences from PLHIV helps us to infer the transmission network5,78. Because HIV-1 has a particularly high evolutionary rate, there is wide genetic diversity within and between hosts. Molecular phylogenetics enables us to quantify pairwise distances between tips on a phylogenetic tree to infer when viruses were likely to have diverged from each other, as in the case of a transmission event. A well-informed molecular HIV phylogeny can provide an estimate of the between host transmission tree4,113. Identifying clusters of closely related HIV sequences has become a common method to detect subpopulations exposed to higher transmission rates than the general population6,75. There are a multitude of methods applied to build clusters from phylogenetic trees and results can vary widely depending on the method8,9. For instance, the chosen genetic distance threshold for cluster detection has a strong effect on which clusters are recovered from the true transmission network76,77. Additionally, cluster detection has been criticized for detecting differences in 45 sampling, rather than differences in transmission in subpopulations81. Furthermore, evaluating whether a patient is within a cluster or not is an over-simplification of the transmission dynamics within clusters. Other attempts have been made to quantify the evolution of pathogen transmission rates across phylogenies using either a Markov-modulated Poisson process82 or a multi-type birth–death branching model process83, however they rely on strong assumptions about the evolution of rates along the tree. Alternatively, the overall transmission history of each tip in a phylogenetic tree can be estimated by its lineage-level diversification rate. Lineage-level diversification rate offers a way of estimating transmission activity without explicitly forming clusters. Viral diversification rate represents the splitting rate of the lineage in the tree, such that a virus will have a high diversification rate if it is evolutionarily similar to many adjacent neighbours94,104. Since transmission of HIV to a new host is equivalent to the formation of a new lineage, diversification rates inferred from phylogenetic trees can serve as estimates for transmission rates generated from patients’ viral genetic sequences. Here, we apply diversification rate for the surveillance of HIV transmission activity and evolutionary dynamics. As a continuous rate, lineage-level diversification rate offers modelling benefits over a logistic regression in which an individual either is or is not a part of a cluster. We hypothesized that higher HIV lineage-level diversification rate of PLHIV in BC would be correlated with known transmission risk factors from cluster analyses, as well as present correlates not previously identified in this cohort. Using HIV sequence data attached to epidemiological, clinical, and demographic data for participants in the Drug Treatment Program, we built a gamma generalized linear model to evaluate significant correlates with HIV transmission in 2018, as estimated using HIV lineage-level diversification rate. A secondary 46 analysis of correlates to the annual change in diversification rate between 2017 and 2018 was also conducted. 3.2 Methods 3.2.1 Study setting and participants Between May 1996 and March 2018, 9,630 participants had drug resistance testing done through their enrolment in the Drug Treatment Program (DTP) at the British Columbia Centre for Excellence in HIV/AIDS (BC-CfE) at St. Paul’s Hospital, Providence Health Care, in Vancouver, BC. The DTP is an ongoing treatment and research cohort42 that provides all medically eligible people living with HIV (PLHIV) in BC with access to HAART and related laboratory monitoring at no direct cost. Eligible are all HIV infected BC residents whether or not they access HAART. Baseline and follow up plasma samples are routinely collected by physician request to evaluate HIV drug resistance and direct optimal use of HAART. Since 2010, the first sample per patient with viral load < 250 copies/mL has been automatically genotyped and did not need a physician request. Notably, not all of the first samples available in this cohort are baseline as migrants from outside the province were also included. Research use samples were also included. Participants in the DTP have been continually recruited since 1992 and follow up continues until death or migration away from BC. Participant enrolment in the DTP occurs at the earliest of either the date of their physician’s first request for medication or at the date of first viral load test. See Supplementary Appendix for details on available patient data. Ethical approval for this study was granted by the University of British Columbia - Providence Health Care Research Ethics Board (H17-01812). No further samples were collected nor were additional patients recruited. 47 3.2.2 Clinical, epidemiological, and demographic variables In addition to sequences, available data for this analysis included sample collection date; date of first viral load; plasma viral load (HIV RNA copies/mL); HIV subtype classifications (generated using the SCUEAL algorithm114); ethnicity; birth year; sex at birth; self-reported risk factors (injection drug user, men who have sex with men, heterosexual contact); having ever tested positive for hepatitis C virus (HCV) infection; having ever had AIDS; if applicable date of mortality and cause of death; and health authority (HA) of patient residence. Regional health authorities in BC are illustrated in Table 3–1; shapefile was obtained from DataBC as a WGS80 projection of NAD83 datum (A.3 Supplementary data and shapefiles). The viral load detection limit for the cross-sectional timepoint in 2018 was 40–10,000,000 copies viral RNA/mL plasma using a Roche COBAS HIV-1 Ampliprep Taqman assay51. HIV sequences and patient information were stored in secure, access-restricted facilities at the BC-CfE in a secure, password protected Oracle database. Patient data were de-identified and doubly anonymized with randomly generated 6-character identifiers. No document was created linking identifiers to the patient. Age was converted to a categorical variable with four levels with similar numbers of PLHIV in each: 41 and under; 42-52; 53-60; and 61 and over. 48 Figure 3-1 Five regional health authorities in British Columbia. 3.2.3 Sequences A total of 36,271 HIV resistance genotype tests from 9,630 PLHIV were completed between May 1996 and March 29, 2018 at the BC-CfE through the DTP. Anywhere from one to 46 plasma tests were run per individual (mean, 3.77; median, 2). The clinical and research staff at the BC-CfE conducted all laboratory work including RNA extraction, PCR, and sequencing. Plasma samples were screened for drug resistance mutations by sequencing the HIV protease and partial reverse transcriptase genes, hereafter referred to as partial pol. The majority of partial pol sequences were 1497 bases (b) in length (n=30,445), encoding protease and the first 400 codons of reverse transcriptase. Another group of sequences (n=5,118) spanned 1017 b, which covered all of protease and the first 240 codons of reverse transcriptase. Finally, a group of 709 sequences were assembled from three partial sequences of 276 b, 297 b, and 363 b, for a total of a 936 b concatenated sequence covering protease and reverse transcriptase codons 24-236, reflecting the progression of HIV testing methodologies. 49 3.2.4 Phylogenetic inference All 36,271 HIV partial pol sequences were aligned to HXB2 reference genome (GenBank Accession #K03455) using MAFFT version 7.310115. The alignment was visualized and curated in AliView V1.17.1116. Insertions and deletions relative to HXB2, as well as amino acids corresponding to WHO recognized drug resistance mutation sites117 were removed from the alignment prior to tree inference. A set of shuffled bootstrap alignments were generated to infer 100 approximate maximum likelihood phylogenetic trees, as implemented in FastTree2.174. Subsequently, trees were pruned to include each patient’s oldest sample and then rooted using root-to-tip regression in ape version 5.0118 using the sampling dates as priors to inform their most recent common ancestor119. 3.2.5 Diversification rate For each tip on each 2018 bootstrap tree, the viral lineage-level diversification rate was calculated and the mean diversification rate across 100 bootstrap trees for each patient at each time interval was carried forward for modelling. In contrast to clade-level measures of diversification98, lineage-level diversification rate is an integrated splitting rate that is unique to each tip in the tree and considers the entire path from root to tip94. It is the reciprocal of the equal splits (ES) metric of evolutionary distinctiveness, which divides the phylogenetic distance of each edge by its daughter edges until the entire tree has been accounted for104. Since ES is an indication of tips’ evolutionary uniqueness, its reciprocal, diversification rate, is an estimate of the extent to which tips are evolutionarily shared. A virus in a rapidly diversifying clade is typified by short branch lengths shared among many tree neighbours, resulting in a low ES value and a high diversification rate value94. Since the viral phylogenetic tree is an estimate of a transmission tree, the diversification 50 rate can be interpreted as an estimate of the between-host transmission rate for epidemiological purposes. The diversification rate for each tip on a rooted bifurcating tree is the reciprocal sum of Ni branch lengths (lj ) from tip i to the root, with each consecutive edge (j) down-weighted by a factor of 1⁄294. 𝐷𝑖𝑣𝑒𝑟𝑠𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝑅𝑎𝑡𝑒2 = 𝐸𝑞𝑢𝑎𝑙 𝑆𝑝𝑙𝑖𝑡𝑠240 = ($ 𝑙&2&4056&/0 )40 The annual change in diversification rate for each participant was calculated as the difference between their 2018 diversification rate (mean of bootstraps) and 2017 diversification rate (mean of bootstraps); newly diagnosed cases between January and March 2018 were assigned a value of zero to emphasize transmission activity among the prevalent population. 3.2.6 Distribution of diversification rates across health authorities and age group As a preliminary investigation, we were interested in whether the differences in diversification rate across age groups were consistent in different Health Authorities. There are five regional health authorities in Vancouver (Figure 3-1). Some participants did not record a health authority because they had unstable housing, were new to the area, or were from out of province. For each health authority separately, we applied a Wilcoxon ranked mean test between each age category with a significance level of 0.05. 3.2.7 Modelling The analysis was restricted to participants who were alive at the time of data generation (sequences from all participants ever in the DTP were used to inform the tree) and had a detectable viral load. Patient epidemiological, clinical, and geographic variables were evaluated in a gamma generalized linear model with a log link. The outcome, lineage-level diversification rate, approximately 51 followed a gamma distribution. Adjusted relative risks were calculated by exponentiating coefficients and confidence intervals. As the goal was to calculate adjusted relative risks for many variables that in some cases were not significantly associated with the outcome, variable selection erred on the side of full inclusivity. 3.3 Results 3.3.1 Patient population There were a total 7,837 PLHIV in BC in March 2018 with available sequences in the DTP dataset who were alive and had a detectable viral load (Table 3–1). Of them, 6,125 were male at birth (81.1% of reported) and 1,351 were female at birth (18.1% of reported). Looking at participants’ age at the time the data was obtained, 1,595 (20.4%) were 41 and under; 2,306 (29.4%) were 42-52; 2,193 (28.0%) were 53-60; and 1,689 (21.6%) were 61 and over. The majority of participants (n=4,172; 53.2%) were located in the Vancouver Coastal Health Authority, followed by Fraser Health Authority (n=1,616; 20.6%), Vancouver Island Health Authority (n=727; 9.3%); unreported health authority (n=672; 8.6%); Interior Health Authority (n=397; 5.1%); and the Northern Health Authority (n=253, 3.2%). Of those who reported to be male at birth, 2,983 (48.1%) were men who have sex with men. As for other risk factors, 2,405 were injection drug users (37.7% of those who reported); 2,561 had been co-infected with HCV previously (36.1% of reported); and 875 had AIDS previously (11.7% of reported). Finally, 3,251 participants self-identified as white (68.3% of reported); 986 as First Nations (20.7% of reported); and 220 as black (4.7% of reported). 52 3.3.2 Distribution of diversification rate by health authority and age category The relationship between age group and lineage-level diversification rate was affected by patients’ health authority of residence (Figure 3-2). Wilcoxon ranked mean tests were used to evaluate whether differences in diversification rate were significantly different between age groups for each health authority. In the Interior Health Authority, there were no significant differences in diversification rate distributions between age categories. In all other health authorities, including unreported, the 41 and under age category exhibited significantly higher diversification rates than the oldest group, 61 and over. Further, in the Northern, Vancouver Island, and Vancouver Coastal Health Authorities, the diversification rates of the 41 and under age category were significantly higher than the 53-60 age group. Only in the Northern Health Authority were diversification rates of the 41 and under group significantly higher than in the 42-52 age group. 53 Table 3–1 The study population included Drug Treatment Program participants in British Columbia with sequence data who were alive on March 29 2018. Characteristic n % of total % of reported All patients 7,837 100.0 - Sex at birth: female 1,351 17.2 18.1 Sex at birth: male 6,125 78.2 81.9 Age category: 61 and over 1,689 21.6 21.7 Age category: 53-60 2,193 28.0 28.2 Age category: 42-52 2,306 29.4 29.6 Age category: 41 and under 1,595 20.4 20.5 Health Authority: Interior 397 5.1 5.5 Health Authority: Vancouver Island 727 9.3 10.1 Health Authority: Fraser 1,616 20.6 22.6 Health Authority: Vancouver Coastal 4,172 53.2 58.2 Health Authority: Northern 253 3.2 3.5 Health Authority: Unreported 672 8.6 - Men who have sex with men: no 2,664 34.0 47.2 Men who have sex with men: yes 2,983 38.1 52.8 Men who have sex with men: unreported 2,105 26.9 - Heterosexual activity: no 5,869 74.9 75.7 Heterosexual activity: yes 1,883 24.0 24.3 Injection drug use: no 3,971 50.7 62.3 Injection drug use: yes 2,405 30.7 37.7 Injection drug use: unreported 1,376 17.6 - Hepatitis C infection ever: no 4,541 57.9 63.9 Hepatitis C infection ever: yes 2,561 32.7 36.1 Hepatitis C infection ever: unreported 374 4.8 - AIDS ever: no 6,601 84.2 88.3 AIDS ever: yes 875 11.2 11.7 Log10(plasma viral load) 7,837 100.0 100.0 Reported ethnicity, white: no 1,511 19.3 31.7 Reported ethnicity, white: yes 3,251 41.5 68.3 Reported ethnicity, First Nations: no 3,776 48.2 79.3 Reported ethnicity, First Nations: yes 986 12.6 20.7 Reported ethnicity, black: no 4,540 57.9 95.3 Reported ethnicity, black: yes 222 2.8 4.7 54 Figure 3-2 The association of age category and lineage-level diversification rate varies by health authority. The difference between mean ranks of diversification rate for different age categories, stratified by health authority, was evaluated using Wilcoxon ranked mean tests, where p > 0.05 (ns), p £ 0.05 (*), p £ 0.01 (**), p £ 0.001 (***), p £ 0.0001 (****). 55 3.3.3 Epidemiological, clinical, and demographic corelates of diversification rate The lineage-level diversification rates in the study population followed a gamma distribution with a right skew and thus a gamma generalized linear model with a log-link equation was applied to estimate adjusted relative risks (Figure 3-3). Participants who were male at birth on average had 1.12 (95% CI, 1.03-1.21) times the adjusted risk of an elevated diversification rate than participants who were females at birth. Those who were 41 years of age or under as well as those who were 42-52 years of age on average had respectively 1.66 (95% CI, 1.52-1.82) and 1.29 (95% CI, 1.19-1.40) times higher diversification rates than the 61 and over age group. Looking at health authorities, participants who lived in the Northern Health Authority had on average 1.31 (95%CI, 1.08-1.59) times higher diversification rate than the control Interior Health Authority group. Those with unreported health authority had on average lower diversification rates than the control group (aRR=0.76, 95%CI, 0.63-0.91). Identifying as MSM was not significantly associated with higher diversification rates in this analysis. However, those who identified as an IDU had significantly higher diversification rates than those who did not (aRR=1.31, 95%CI, 1.19-1.45); similarly, those who ever had HCV infection had on average 1.23 (95%CI, 1.13-1.34) times higher diversification rates than those who had not, adjusted for other risk factors. Individuals who ever had AIDS on average had lower diversification rates (aRR=0.76, 95%CI, 0.70-0.83). The most recent plasma viral load was also significantly associated with diversification rate – for every 1-unit increase in the log10(viral load), individuals had on average 1.14 (95%CI, 1.11-1.17) times higher diversification rates, adjusted for other factors. Participants who identified as black had on average lower diversification rates (aRR=0.62, 95%CI, 0.52-0.75). 56 Figure 3-3 A summary of the adjusted relative risks for patient attributes associated with HIV lineage-level diversification rates. Risks were estimated using a multiple gamma generalized linear regression model with a log-linker equation. A parallel analysis was conducted using the annual change in diversification rate between 2018 and 2017 as the outcome (Figure 3-4). Similar to the analysis of diversification rate, the 41 and under age group on average had significantly higher changes in diversification rate than the 62+ age group (aRR=1.90, 95%CI, 1.52-2.38), as did the 42-52 age group (aRR=1.45, 95%CI, 1.19-1.77). Interestingly, of all the health authorities, participants in Vancouver Island had the highest risk of large changes in diversification rate relative to the Interior Health reference group (aRR=1.50, 95%CI, 1.02-2.16). Living in the Northern Health Authority, being an IDU, ever having HCV, ever having AIDS, viral load, and black ethnicity were no longer significant. A new 57 predictor arose, however – those who identified as white ethnicity on average had 1.45 (95%CI, 1.13-1.84) times higher changes in diversification rate than those who did not identify as white. Figure 3-4 A summary of the adjusted relative risks for patient attributes associated with annual changes in HIV lineage-level diversification rate between 2018 and 2017. Risks were estimated using a multiple gamma generalized linear regression model with a log-linker equation. 3.4 Discussion In this cross-sectional analysis of PLHIV in BC in 2018, we have demonstrated that high HIV lineage-level diversification rates were significantly associated with being male, young, living in the Northern Health Authority, being an IDU, ever having been infected with HCV, or having a high most recent viral load, while having a low diversification rate was associated with ever having AIDS, self-reporting black ethnicity, or not reporting a health authority. In a secondary analysis, having a high annual change in diversification rate was associated with being young, living in the 58 Vancouver Island Health Authority, and being white. Diversification rate could be used to complement existing phylogenetic public health programs to prioritize enhanced care and connection to treatment for those who may be at an elevated risk of participating in ongoing transmission. Participants living in the Northern Health Authority had significantly higher diversification rates than the reference group, despite having the smallest population of PLHIV within the sequenced DTP cohort. It is unclear whether this is due to this population having a harder time accessing health care resources for testing, treatment, and follow-up, or whether this could indicate more complete sampling coverage in that area81. By contrast, participants living in the Vancouver Island Health Authority had significantly higher annual changes in diversification rate, which could be more indicative of most recent transmission activity among those sampled. Being young (particularly under 41 or 42-52) was a risk factor for HIV transmission in every health authority except Interior Health, also when considering the annual change in diversification rate. This could indicate that these groups have a harder time connecting with care. Conceptually, it makes sense that younger people would be more active socially, sexually, and may be more likely to engage in risk-taking behaviour34. Although it was demonstrated that the effect of age category on diversification rate was modified by health authority, no effect modification term was included in the model to err on the side of interpretability of adjusted relative risks as that would have included another 24 coefficients into the model (six health authorities and 4 age categories). The age groups were made on the basis of having four groups with similar numbers of PLHIV in each and results may have differed slightly had these arbitrary cut-offs been different or had there been more than four categories. 59 Although being MSM was not deemed to be a significant risk factor for HIV transmission in our model and was a protective factor in the cluster analysis, this could have been due to lower than average sampling coverage of PLHIV who are MSM. In Vancouver Coastal Health Authority, in particular, the majority of new HIV cases in recent years have arisen in the West End of Vancouver among young MSM31. It is possible that a high percent of the undiagnosed HIV+ population are young MSM. A low sampling coverage from this population would partially explain why both this analysis and the cluster analysis by Poon et al. did not identify that being MSM put someone at a significantly higher risk of HIV transmission81. That IDU and participants ever infected with HCV (which are colinear in the first place) were significant when evaluating diversification rate, but not when evaluating the annual change in diversification rate was consistent with what we know about the HIV epidemic in Vancouver. While in the 1990’s, injection drug use was a primary risk factor for HIV and HCV transmission28, harm reduction measures targeted at preventing HIV transmission among IDU such as InSite, North America’s first safe injection site32, have markedly reduced new HIV cases in this population31. We compared the significant associations identified by our diversification rate analysis to those obtained in the Poon et al. cluster analysis using an earlier version of the dataset22. It is relevant to note that there were several variables in either analysis that were not evaluated in the other (i.e. drug resistance history and transmitted drug resistance in the cluster analysis; and ethnicity groups and age category at time of data acquisition in the diversification rate analysis). Our methods also differed in that we excluded those who had died and that we calculated adjusted 60 relative risks while Poon et al. calculated univariate odds ratios. Despite these differences, it is valid to qualitatively compare these results as the trees and participant sampling distributions were nearly identical. Both our diversification rate analysis and the cluster analysis found that injection drug use and ever having had HCV were significantly associated with high transmission (aRR>1; OR>1), while ever having had AIDS was associated with lower transmission on average (aRR<1; aRR<1). Our analysis additionally identified being male, young, having a high most recent viral load, and living in the Northern Health Authority as risk factors for elevated transmission (aRR>1), as well as finding that participants who identified as black and who did not report a health authority generally had lower transmission activity (aRR<1). By contrast, the cluster analysis identified that females were more likely to be in clusters than males (OR>1). They also found that being older at baseline, having a drug resistance history, and being MSM were all associated with lower probability of clustering (OR<1), while having a high baseline viral load, showing transmitted drug resistance, and heterosexual activity were associated with higher probability of clustering (OR>1). In the future we will repeat this comparison with an identical dataset, phylogeny, participant restrictions, and modelling strategy to emphasize differences in methodologies alone, as well as assess the potential of combining clustering and diversification rate. Another notable limitation of our analysis is that we disregard those who are unsampled in the population. Our conclusions are thus only fully generalizable to those who are connected with care and were represented in this study. A further limitation is that lineage-level diversification rate may be partially subject to a similar challenge as clustering in that differences between sub-populations may reflect sampling differences more readily than transmission differences81. In future work, we will be comparing the sensitivities of diversification rate and clustering in response 61 to changes in sampling coverage. A final limitation is that the annual change in diversification rate between 2018 and 2017 only included samples until March 29 2018. For clarification, an individual can have a high diversification rate without themselves being an active infector, by virtue of their virus being phylogenetically proximate to the virus of someone who is actively transmitting HIV. A group risk does not necessarily equate to an individual risk and diversification is a shared risk within groups of proximate lineages. Lineage-level diversification rate should not be used to support any trial evidence that criminalizes HIV non-disclosure. Phylogenetics are inherently uncertain and should be interpreted thus. It is imperative that we translate phylogenetic findings into real changes in the way public health care is provided. We recommend the groups we have identified as being at risk have additional outreach and services be made immediately available. For instance, normalize and expand point-of-care testing for young people in all health authorities; engage PLHIV with their care to keep them adherent with treatment; and for PLHIV within risk groups, support their partners to engage with PrEP. 3.5 Conclusion Studying epidemiological risk factors for HIV transmission using lineage-level diversification rate permitted us to corroborate known risk factors and illuminate novel associations, demonstrating its potential applicability in the field of infectious disease epidemiology. 62 Chapter 4: Viral diversification rate identifies HIV transmission hot spots Angela McLaughlin1,2, Ph.D., Paul Sereda1, Rolando Barrios1, Ph.D., Chanson J. Brumme1, Ph.D., Zabrina L. Brumme1,3, Ph.D., Julio S.G. Montaner1,4, M.D., Natalia Oliveira1, and Jeffrey B. Joy1,4, Ph.D. 1BC Centre for Excellence in HIV/AIDS, St. Paul’s Hospital, Vancouver, BC 2School of Population and Public Health, University of British Columbia Vancouver, BC 3Faculty of Health Sciences, Simon Fraser University, Burnaby, BC 4Division of Infectious Diseases, Department of Medicine, University of British Columbia, Vancouver, BC Identifying key populations at high risk of ongoing HIV transmission is critical for prioritizing limited health care resources, preventing new HIV infections, and achieving the UNAIDS 90-90-90 Target. Lineage-level diversification rates inferred from viral phylogenetic trees provide granular estimates of HIV transmission; each tip’s lineage-level diversification rate is its weighted splitting rate from tip to root. To identify HIV-1 transmission foci in British Columbia, Canada, we inferred viral diversification rates from phylogenetic trees of 36,271 HIV-1 sequences from 9,630 anonymized patients. Diversification rates were combined with patients’ epidemiological and clinical data and aggregated by area of patients’ residence to predict the spatial distribution of new HIV cases across British Columbia between 2008 and 2018. The predictive power of the model was compared with a phylogenetically uninformed model. For every additional participant in a census tract with an annual change in diversification rate greater than one, there were an average of 1.10 (95% CI, 1.07-1.13) times more new HIV cases in the subsequent year. Further, for every 1-unit increase in the mean of the top five diversification rates, the number of new HIV cases increased by on average 1.38-fold (95% CI, 1.28-1.49). In a blind prediction of 2018 cases, viral diversification rate metrics improved the model’s specificity (0.46; 0.41), accuracy (0.50; 0.46), and correlation of predicted and observed values (0.42; 0.16). Geographically-aggregated viral diversification rates were predictive of new HIV cases in small geographic areas and they improved the fit of a predictive model relative to epidemiological and clinical correlates alone. 63 4.1 Introduction In an era of limited public health resources, the efficient and equitable allocation of health services is increasingly important. Specifically, identifying geographic areas that are at high risk of ongoing transmission of HIV is critical for the efficient allocation of prevention and treatment services to meet the UNAIDS 90-90-90 Target35. Although highly active antiretroviral therapy (HAART) has led to sustained decreases in morbidity, mortality, and transmission for people living with HIV (PLHIV)30, foci of high transmission remain, even in developed countries6,22. Secondary analyses of pre-existing HIV resistance genotyping datasets collected routinely to personalize patients’ ART enable us to study viral transmission as it relates to patient characteristics and epidemiological outcomes22. For rapidly evolving pathogens like HIV, a well-informed molecular phylogeny can provide an estimate of the between-host transmission tree4,113. Phylogenetic clustering of HIV sequences is commonly used to detect subpopulations exposed to higher transmission rates in order to prioritize services such as pre-exposure prophylaxis (PrEP) and enhanced partner contact tracing6,7,75. However, no standard cluster detection method has been agreed upon and results can vary dramatically81,112. Alternatively, the overall transmission history of each tip in a phylogenetic tree can be estimated by its lineage-level diversification rate, which is a tip’s splitting rate across the entire tree path, weighted from tip to root94,104. Since transmission of HIV to a new host is equivalent to the formation of a new lineage, diversification rates inferred from phylogenetic trees can serve as a proxy for transmission rates inferred from patients’ viral genetic sequences. In previous work, we have highlighted that people who inject drugs (PWID) and men who have sex with men (MSM) 64 exhibited significantly higher HIV diversification rates and that treatment experienced lineages displayed dramatically reduced HIV diversification rates, providing an independent validation of Treatment as Prevention® (TasP®) strategies95. Aggregating viral diversification rates by patients’ geography of residence, therefore, should illuminate areas with concentrated transmission activity. Using longitudinal data collected from participants in the Drug Treatment Program (DTP) through the BC Centre for Excellence in HIV/AIDS (BC-CfE), we sought to model the spatiotemporal distribution of HIV lineage-level diversification rate in BC in order to predict where new HIV cases arise over time. We hypothesized that geographically-aggregated viral diversification rate would be a significant predictor of the number of new HIV cases, lending significantly higher predictive power than a model that considers only clinical and epidemiological patient characteristics. 4.2 Methods 4.2.1 Study setting and participants Between May 1996 and March 2018, 9,630 participants had drug resistance testing done through their enrolment in the Drug Treatment Program (DTP) at the BC-CfE at St. Paul’s Hospital, Providence Health Care, in Vancouver, BC. The DTP is an ongoing treatment and research cohort42 that provides all medically eligible people living with HIV (PLHIV) in BC with access to HAART and related laboratory monitoring at no direct cost. Eligible are all HIV infected BC residents whether or not they access HAART. Baseline and follow up plasma samples are routinely collected in all participants to perform HIV genotyping to evaluate HIV drug resistance to direct optimal use of HAART. Participants in the DTP have been continually recruited since 1992 and follow up continues until death or migration away from BC. Participant enrolment in the DTP 65 occurs at the time of their physician’s first request for medication. See Supplementary Appendix for details on available patient data. Ethical approval for this study was granted by the University of British Columbia - Providence Health Care Research Ethics Board (H17-01812). No further samples were collected nor were additional patients recruited. 4.2.2 Phylogenetic inference A total of 36,271 HIV resistance genotype tests from 9,630 PLHIV were completed between May 1996 and March 29 2018 at the BC-CfE through the DTP. Anywhere from one to 46 plasma tests were run per individual (mean, 3.77; median, 2). Plasma samples were screened for drug resistance mutations by sequencing the HIV protease and partial reverse transcriptase genes, hereafter referred to as partial pol. See Supplementary Appendix for details on the sequence data. Sequences were aligned to the HXB2 reference genome (GenBank Accession #K03455) using MAFFT version 7.310115. The alignment was visualized and curated in AliView V1.17.1116. Insertions and deletions relative to HXB2, as well as amino acids corresponding to WHO recognized drug resistance mutation sites117, were removed from the alignment prior to tree inference. A set of shuffled bootstrap alignments were generated to infer 100 approximate maximum likelihood phylogenetic trees (in units of substitutions per site), as implemented in FastTree2.174. Subsequently, trees were pruned to include each patient’s oldest sample and then rooted using root-to-tip regression in the R package ape version 5.0118 using the sampling dates to inform their ancestry119. To evaluate temporal changes, rooted trees were then pruned to generate 100 trees per time interval (Figure A-3). Each time interval tree included patients’ oldest sequences sampled that would have been available in that year, such that each year the trees only had tips added for newly infected (or newly connected with care) participants. Alternative tree-building algorithms and respective diversification rate calculations were explored in Appendix A. 66 4.2.3 Diversification rate For each tip on a rooted bifurcating (time interval bootstrap) tree, the viral lineage-level diversification rate (DRi) was calculated as the reciprocal sum of Ni branch lengths (lj ) from tip i to the root, with each consecutive edge (j) down-weighted by a factor of 1⁄294. 𝐷𝑖𝑣𝑒𝑟𝑠𝑖𝑓𝑖𝑐𝑎𝑡𝑖𝑜𝑛 𝑅𝑎𝑡𝑒2 = 𝐸𝑞𝑢𝑎𝑙 𝑆𝑝𝑙𝑖𝑡𝑠240 = ($ 𝑙&2&4056&/0 )40 The annual change in diversification rate for each lineage over time was calculated as the difference between their mean diversification rate across bootstraps for a given year and the preceding year, with a value of zero for the first year a patient entered the cohort. A change in diversification rate indicates that a new patient’s virus has been added to the tree in close phylogenetic proximity to a viral tip. We described the diversification rate in more detail, characterized its distribution over time, and evaluated the robustness of our estimates to different tree-building algorithms in Appendix A. 4.2.4 Predictive modelling Patient attributes and lineage-level diversification rates were aggregated by their census tract of residence. A census tract is a small geographic area with generally under 10,000 people within a census agglomeration with a population of 50,000 or more in the previous census120. See Appendix A for further details on patient attributes, aggregated variables, and geographic areas of analysis. The primary outcome of this study was defined as the number new HIV cases in each census tract in the subsequent year, as estimated by date of first detectable viral load (see Appendix A for outcome validity discussion). Since this is a discrete count with a large percentage of zero values and some over-dispersion (Figure A-2), a zero-inflated negative binomial (ZINB) model was 67 deemed most appropriate121. A simple supervised machine learning algorithm was applied to train and test the model (workflow in Figure A-7). The spatiotemporal distribution of the outcome and predictors were visualized (downtown Vancouver, Figure 4-2). To assess how informative diversification rate measures were in predicting new HIV cases, the full ZINB model was reduced to a nested model excluding any phylogenetic measures (referred to as reduced ZINB). The goodness-of-fit of these models, a full zero-inflated Poisson (ZIP), and a full Hurdle model were compared in Table A-3. The predictive fit of the full ZINB and reduced ZINB were compared in terms of their Pearson’s correlation coefficient for observed values greater than zero; sensitivity; specificity; positive predictive value; negative predictive value; and accuracy (Figure 4-3; Figure A-6). Finally, the model was used to blindly predict priority areas where HIV transmission was highest in 2018 based on 2017 predictors. 4.3 Results 4.3.1 Study population Of the 14,235 cumulative DTP participants by June 30 2018, 9,630 individuals had sequence data by March 29 2018, and 6,944 of them reported their baseline census tract of residence or postal code and were included in the final analysis (Table 4–1). The restricted population for the analysis was comparable to the DTP population with sequences in most regards including the distribution of sexes (male: 81.3% among DTP participants with sequences, 83.3% among restricted population); dominant HIV subtype (subtype B: 92.4%, 92.5%); median baseline age (38, 38); and median baseline diversification rate (52.5 (subs/site)-1, 52.0 (subs/site)-1). The study population slightly underrepresented participants from the Vancouver Island (10.1%, 7.2%), Interior (5.7%, 2.9%), and Northern (3.6%, 1.9%) Health Authorities, and First Nations people (22.2%, 18.9%). 68 Additionally, the median baseline viral load was somewhat lower in the DTP population with sequences than among those who reported a census tract (104.5 copies/mL, 105.3 copies/mL). Table 4–1 Characteristics of the Drug Treatment Program cohort with sequences compared to the study population restricted to those who have reported their census tract of residence. DTP cohort with sequences Study population Parameter Characteristic N* Total % N* Total % Total participants 9,630 9,630 100 6,944 6,944 100 Sex Male 8,861 7,208 81.3 6,634 5,525 83.3 Female 8,861 1,592 18.0 6,634 1,053 15.9 Male-Female 8,861 51 0.6 6,634 47 0.7 Female-Male 8,861 10 0.1 6,634 9 0.1 Sex at birth Male 8,861 7,259 81.9 6,634 5,572 84.0 Female 8,861 1,602 18.1 6,634 1,062 16.0 Subtype B 9,630 8,895 92.4 6,944 6,420 92.5 C 9,630 361 3.7 6,944 254 3.7 AE 9,630 163 1.7 6,944 118 1.7 A 9,630 85 0.9 6,944 59 0.8 AG 9,630 60 0.6 6,944 43 0.6 Other 9,630 34 0.4 6,944 26 0.4 D 9,630 32 0.3 6,944 24 0.3 Health Authority Vancouver Coastal 8,459 4,923 58.2 6,589 4,217 64.0 Fraser 8,459 1,894 22.4 6,589 1,581 24.0 Vancouver Island 8,459 858 10.1 6,589 475 7.2 Interior 8,459 478 5.7 6,589 191 2.9 Northern 8,459 306 3.6 6,589 125 1.9 Self-reported risk factors Injection drug use 5,500 3,175 57.7 3,948 2,146 54.4 Men who have sex with men 6,874 3,487 50.7 5,051 2,783 55.1 Heterosexual contact 6,874 2,201 32.0 5,051 1,562 30.9 Other 6,874 308 4.5 5,051 214 4.2 Exposure to blood products 6,874 225 3.3 5,051 161 3.2 Self-reported ethnicity White 5,753 3,944 68.6 4,284 2,998 70.0 First Nations 5,753 1,275 22.2 4,284 811 18.9 Asian 5,753 435 7.6 4,284 378 8.8 Black 5,753 247 4.3 4,284 192 4.5 Hispanic 5,753 221 3.8 4,284 176 4.1 Clinical illness Ever had hepatitis C virus 8,318 3,268 39.3 6,247 2,265 36.3 Ever had AIDS 8,861 1,074 12.1 6,634 824 12.4 Died 7,231 1,527 21.1 5,283 966 18.3 Baseline measures N* Median IQR N* Median IQR Baseline diversification rate (subs/site)-1 9,630 52.5 34.4 - 89.1 6,944 52 34.3 - 88.3 Baseline log10 (viral load) log10(viral RNA/mL plasma) 9,630 4.5 3.8 - 5.0 6,944 5.3 4.5 - 7.0 Baseline age 9,624 38 31 - 46 6,916 38 31 - 45 *N = the number of participants who provided information for a given characteristic 69 4.3.2 Diversification rate The distribution of lineage-level diversification rates across a representative approximate maximum likelihood tree from 2017 is presented in Figure 4-1. Lineages with rapid and extensive branching are typified by high diversification rates. We showed that our diversification rate calculations were robust to tree-building methods and genomic region (Table A-1). Although there was a slight increase in both mean and median diversification rate over time attributable to sampling from a larger tree (Figure A-4), a Kruskal-Wallis test revealed that the population-level distribution of diversification rates did not significantly change across the study period (p=0.385). Figure 4-1 A representative bootstrap approximate maximum likelihood phylogenetic tree with branches coloured by the lineage-level diversification rate. Tree includes baseline sequences available for all participants as of March 29 2018. 70 4.3.3 Predictive modeling The unadjusted relationships of aggregated variables with the outcome were calculated to inform their model inclusion, as summarized in Table A-2. Among the highest Pearson correlation coefficients were new HIV cases in the previous year (r=0.74, p<0.001), the total number of PLHIV with changes in diversification rate ³ 1 (subs/site)-1 (r=0.72, p<0.001), the sum of PLHIV for whom all risk info was known (r=0.71, p<0.001), the sum of PLHIV (r=0.69, p<0.001), and the mean of the top 5 diversification rates (r=0.48, p<0.001). After identifying a final ZINB model, a likelihood ratio test confirmed that it fit the data significantly better than its Poisson equivalent (p<0.001), Hurdle equivalent (p<0.001), and the reduced ZINB without phylogenetic measures (p<0.001) (Table A-3). Further, the AIC, BIC, and log-likelihood values for the aforementioned models were all lower for the full ZINB model. The odds ratios and relative risks associated with significant predictors in the final model show that the mean of the top 5 ln(diversification rates) was significant in both parts of the model (Table 4–2). The adjusted odds of a census tract having greater than zero new HIV cases in the subsequent year was 3.10 (95% CI, 1.61-5.97) times higher for every 1-unit increase in the mean of the top 5 ln(diversification rates). Given that a census tract had greater than zero new HIV cases, for every 1-unit increase in the mean of the top 5 ln(diversification rate) of a census tract, the adjusted risk of new HIV cases in the subsequent year was 1.38 (95% CI, 1.28-1.49) times higher. It was further of interest that for every additional PLHIV with change in diversification rate ³ 1 (subs/site)-1 in a census tract, the adjusted risk of new HIV cases in the subsequent year was 1.10 (95% CI, 1.07-1.13) times higher. The spatiotemporal distribution of new HIV cases in downtown 71 Vancouver across the study period was heterogenous and mirrored trends among the significant predictor variables (Figure 4-2). Table 4–2 The final zero-inflated negative binomial (ZINB) predictive model. Within the two-part model, the binomial model predicts whether there were greater than zero new cases in any census tract and the negative-binomial model predicts the number of cases, if cases were greater than zero. Binomial {0,>0} model Adjusted odds ratio 95% CI p-value Total PLHIV 0.81 0.74 - 0.88 <0.001 Mean of top 5 ln(diversification rate) 3.10 1.61 - 5.97 <0.001 Negative binomial {count if >0} model Adjusted relative risk 95% CI p-value Total new cases 1.07 1.03 - 1.11 <0.001 Mean of top five ln(diversification rate) 1.38 1.28 - 1.49 <0.001 Total # PLHIV with change in diversification rate ³ 1 1.10 1.07 - 1.13 <0.001 Comparing the predictive fit of the ZINB model to its reduced equivalent for the testing dataset, the mean Pearson correlation coefficient of predicted and observed values for the full ZINB model of 0.77 was higher than that of the reduced ZINB model at 0.66 (Figure 4-3; Table A-4). Interestingly, the reduced ZINB model had a higher mean sensitivity than the full ZINB model (0.75, 0.66), but the full ZINB model had a higher mean specificity (0.68, 0.61) due to the reduced ZINB producing more false positives. The mean accuracy was the same for both models (0.67). Of the 201 new HIV cases with sequences in 2018, only 64 had reported census tracts or postal codes in BC by February 4 2019, a markedly lower percent than previous years (Table A-5). Regardless, observed new HIV cases in 2018 were compared to blind predictions from the model using 2017 predictors (Figure 4-3; Figure A-6). The Pearson correlation between observed and predicted number of HIV cases for the full ZINB model was markedly better than the reduced model (0.42, 0.16). The sensitivity for both models was quite high, although the reduced model 72 had a higher sensitivity of 0.90 compared to 0.79 in the full model. The full ZINB model had a higher specificity than the reduced model (0.46, 0.41) and a higher accuracy than the reduced ZINB model (0.50, 0.46). Figure 4-2 Heat maps of the epicenter of the HIV epidemic in BC, downtown Vancouver, on odd study years between 2008 and 2017. Even years were also included as time points in the study but were excluded here for conciseness. The outcome, (a) total new HIV cases, was predicted by (b) total prevalent cases of PLHIV, (c) sum of PLHIV with annual changes in diversification rate (DR) ³1 (subs/site)-1, and (d) mean of top 5 ln(diversification rates, DR). 0km 2km 4km2009 2011 2013 2015 2017TotalNew Cases051015200km 2km 4km2009 2011 2013 2015 2017Sum PLHIVwith DRchange >= 1051015200km 2km 4km2009 2011 2013 2015 2017Mean (top5 ln(DR))45670km 2km 4km2009 2011 2013 2015 2017SumPLHIV0100200300400(a) (b) (c) (d) 73 Figure 4-3 A comparison of the predictive power of the full ZINB and ZINB without phylogenetic measures (reduced ZINB). Boxplots summarize median and interquartile range values for three separate testing years (2008, 2011, and 2014), while the diamonds illustrate the predictive values for the blind prediction of 2018 new HIV cases. Criteria considered for predictive power include the Pearson’s correlation coefficient (for observed cases > 0 to remove the effect of zero-inflation); sensitivity; specificity; positive predictive value; negative predictive value; and accuracy. 4.4 Discussion Our findings revealed that geographically-aggregated HIV lineage-level diversification rate, supplemented with clinical and epidemiological data, better predicted the number and location of future new HIV cases across BC than clinical and epidemiological data alone. Such predictive models of where new HIV cases will arise in the future are valuable to public health authorities for prioritizing the provision of treatment and prevention services to areas of greatest need. Several measures of aggregated diversification rate were significantly predictive of where new HIV cases arose, even after correction for HIV seroprevalence and incidence in the previous year. This supports the utility of lineage-level diversification rates as estimates of HIV transmission and suggests that the communities where actively transmitting PLHIV reside tend to 74 be those where new HIV cases arise. Though initially counterintuitive, increases in the prevalent PLHIV population of a census tract resulted in a decreased risk of new HIV cases in the subsequent year. This was likely because the census tracts with the highest prevalence are well known to public health officials and already have significant prevention and treatment programs in place. Speaking to the predictive power of the models, the full model exhibited an improved correlation of observed and predicted new HIV cases and a higher specificity in both the testing and blind prediction data relative to the reduced model, whereas the reduced model had a higher sensitivity than the full model. The full model tended to produce fewer false positives and more false negatives than the reduced model. While the differences in the models’ predictive powers were modest, they demonstrated that the inclusion of diversification rate measures was informative to where and how many new HIV cases arose in the subsequent year. Further optimization of variable selection could improve its predictive power further. The 2018 blind predictions were somewhat hindered by the low percent of new HIV cases with sequences who also had census tract information; following up with participants to assemble this data would be worthwhile. Although restricting analyses to individuals who resided in a census tract potentially could have limited the generalizability of the results due to sampling bias against those more likely to have a reported residence, we found that participants with reported census tracts were comparable to the DTP group with sequences available. A greater sampling bias was likely derived from those who were disconnected from care entirely. It was recently estimated that in 2014, 83% of PLHIV in BC were diagnosed35. The extent that undiagnosed PLHIV contribute to new infections is unknown, however it is possible that this fraction might disproportionately contribute to HIV transmission as a result of unknown HIV status, lack of viral suppression, and potentially 75 participating in high risk behaviors123. Efforts should be made to engage hard-to-reach individuals in care in a way that is non-threatening, inclusive, and feasible. Other attempts have been made to quantify the evolution of pathogen transmission rates across phylogenies82,83. These methods fundamentally differ from the methodology described here as these tree-wide estimates assume that the transmission rates evolve as a discrete state according to either a Markov-modulated Poisson process82 or a multi-type birth–death branching model process83, in contrast to our lineage-level estimates, which do not rely on assumptions about the evolution of rates. It is imperative that researchers in genomic epidemiology protect patients’ privacy and confidentiality by mitigating the risk of patient identification89. Efforts were made in this study to protect patients’ privacy and confidentiality by merging census tracts until no fewer than three PLHIV resided in the final merged census tract at any point in time. Another reason for aggregating multiple census tracts comes from the fact that any rates computed using small populations at risk are relatively unstable124. In future work, we intend to validate the model on other HIV-afflicted geographies with lower sampling coverage. Simulating an epidemic could also be useful to assess the effect of random down-sampling on the distribution and interpretation of diversification rates in comparison to cluster detection methods. Further, calculating the lineage-level diversification rates of individuals within clusters could increase the resolution of existing prioritization schemes. Lineage-level diversification rate should prove to be a useful metric for other rapidly-evolving pathogens as well. 76 4.5 Conclusions Viral lineage-level diversification rate is a useful estimate of HIV transmission and when aggregated by patients’ geography of residence revealed high priority transmission hot spots. Phylogenetic analyses of pre-existing HIV genotyping data could refine predictions of where public health resources will have the greatest impact. 77 Chapter 5: Conclusion 5.1 Overall analysis and integration of the research and conclusions Together, my analyses of patient-level and geographically-aggregated viral diversification rates among PLHIV informed which patient attributes were associated with high viral transmission rates and where HIV cases were likely to occur in BC. In the cross-sectional analysis of epidemiological correlates, having a high HIV diversification rate was positively associated with being male, younger, using injection drugs, ever having hepatitis C virus, having a high recent viral load, and residing within the Northern BC Health authority (p<0.001). In contrast, having ever had AIDS and identifying as black were both significantly associated with lower diversification rates (p<0.001). A secondary analysis suggested that higher annual changes in diversification rate were associated with being younger, white, and residing on Vancouver Island (p<0.001). However, the dataset used in the analysis only had sequences up until March 29 2018. The geographic model indicated that HIV diversification rates in conjunction with epidemiological measures were predictive of new HIV cases across census tracts in BC in the subsequent year. This model had an improved specificity, accuracy, and correlation of predicted and observed values compared to the reduced model lacking phylogenetic information for both the testing and blind prediction datasets. Improvements in model fit with phylogenetic metrics were modest, but the diversification rate did add relevant signal to the model. Likely an alternative variable selection algorithm and machine learning method could have been optimized further, however this serves as a preliminary step in that direction. Outside the scope of HIV transmission analyses, predictive geographic models for health outcomes in unique communities have wide scope for the equitable allocation of health care resources. 78 5.2 Re-visiting original thesis hypotheses • HIV lineage-level diversification rate of PLHIV in BC is significantly associated with known transmission risk factors such as injection drug use and being young, in addition to previously unknown putative risk factors. o We failed to reject the null hypothesis that being male, young, having a history of injection drug use, ever having had HCV, having a high most recent viral load, or living in the Northern health authority were unrelated to the risk of having higher than average lineage-level diversification rate. o Likewise, those who ever had AIDS, did not report a health authority, or identified as black had significantly lower lineage-level diversification rates than the average PLHIV population. • Geographically-aggregated HIV lineage-level diversification rate is a significant predictor of where new HIV cases arise in the future. o Multiple aggregated measures of lineage-level diversification rate were correlated with new HIV cases in the subsequent year. o The inclusion of diversification rate measures in a predictive model of new HIV cases in merged census tracts between 2008 and 2017 led to a significantly higher specificity, accuracy, and correlation of predicted and observed values in the testing dataset relative to a reduced model containing only epidemiological predictors. • The predictive model trained and tested using the 2008-2017 data is predictive of ongoing data patterns through 2018. 79 o The blind predictions of new HIV cases in 2018 using the full ZINB model demonstrated a higher specificity (46%, 41%), accuracy (50%, 46%), and correlation (r=0.42, r=0.16) than the reduced ZINB, although the reduced ZINB had a higher sensitivity than the full ZINB (90%, 79%). 5.3 Overall significance and contribution These findings represent significant contributions in the field, as the relationship between diversification rate and new HIV cases had yet to be evaluated. We showed that incorporating phylogenetic measures into a geographic model of new HIV cases in the subsequent year added predictive power to a model excluding phylogenetic measures. This suggested that lineage-level diversification rate is a valid proxy for transmission activity and is an independent predictor of where and how many new HIV cases occur. While the cross-sectional analysis of diversification rate in relation to epidemiological, clinical, and demographic variables yielded similar results to a previous cluster analysis using an earlier version of the dataset, there were several novel associations identified. Namely, those who were currently young, male, had a high most recent viral load, or resided in the Northern health authority had a higher risk of elevated diversification rates, while those who identified as black had on average lower diversification rates. Our secondary analysis of annual changes in diversification rate between 2018 and 2017, a first of its kind, suggested that those who lived on Vancouver Island, were young, or white had significantly higher risks of elevated diversification rates. This analysis was done with the caveat of only having viral sequences up to March 2018 and thus was only representative of the dynamics of the local epidemic in the first quarter of the year. 80 5.4 Strengths and limitations The primary strengths of this research were its translatability and novelty. The findings of both Chapter 3 and Chapter 4 included recommendations for enhanced treatment and prevention services for groups and geographies at risk of ongoing transmission. Additionally, the heat maps of significant predictors and the outcome are easily interpretable for public health officials. As for novelty, these analyses were a unique combination of phylogenetics, epidemiology, geography, and modelling that leaves room for improvement. Our methodologies to apply lineage-level diversification rate to study the evolutionary dynamics of HIV set a new precedent that could be applied to other rapidly evolving pathogens, as well. One of the greatest limitations of the studies was the inherent sampling bias at successive levels, including only those in the Drug Treatment Program, those with sequences, and then for the geographic analysis, those with a reported census tract. At each level, we excluded individuals who may be systematically different from those who we have included and thus our results are truly only generalizable to the population in the DTP with sequences and an address. Another limitation of our analyses is that they assume the partial pol gene used to calculate diversification rates is reflective of evolution across the entire HIV genome. The correlations to diversification rates calculated using integrase were only possible for a small subset of participants and improved with sampling coverage, but could have been higher. Of course, there is the large assumption that the approximate maximum likelihood trees used to inform the diversification rates are somewhat close estimations of the transmission network. For the spatial analysis, we aggregated patient attributes by their census tract of residence at baseline, leading to two big assumptions – first, that patients continue living in the same census tract and second, that people transmit virus in the same area they live in. This affects the interpretation of the conclusions. An additional limitation in the 81 spatial analysis was the use of date of first detectable viral load as an estimate of the date of infection. Of course, there is a delay between infection and diagnosis and a further delay between diagnosis and measuring viral load. However, because new HIV cases were evaluated on an annual basis, only patients with delays in diagnosis and viral load measurement spanning multiple years would have been misrepresented as a new case in the wrong year. New HIV cases also included first detectable viral load from new migrants to BC who may have been previously infected and may be treatment experienced, further limiting our interpretation of new HIV cases as a proxy for new HIV infections. A final notable limitation of these analyses was that the predictive model that included diversification rate metrics had only modest improvements in its specificity, correlation, and accuracy relative to its reduced form, which had a higher sensitivity than the full model. This suggested that while diversification rate did lend additional useful signal to the model, there is still room for improvement in variable selection and coefficient estimation if this model were to be applied in a public health setting. 5.5 Potential applications of the research findings Knowing in what groups and where new HIV cases will likely arise in the future will help us to allocate resources to the people who need them most. Based on our cross-sectional analysis, we believe a greater effort should be put forth in all health authorities to engage young people with prevention (PrEP, education, condoms, partner tracing and contact) and treatment (HAART, adherence follow-up), as well as normalizing that discussion about communicable disease risks with young and old groups alike. We further recommend that the areas with the highest predicted number of new HIV cases in the West End of Vancouver receive consideration as sites for enhanced harm reduction programs such as PreP expansion, outreach education campaigns, condom distribution, and perhaps additional safe injection sites. At the very least, this serves as an 82 additional piece of information for public health authorities to use to inform their decision-making and policy. This is just one example of a predictive geographic model that incorporates epidemiological, clinical, and transmission to predict a health outcome. The potential applications to more accurately map where health care resources are needed most are vast. Further method development to partner lineage-level diversification rate with existing HIV cluster detection methods will provide further resolution to prioritized care strategies and corroborate cluster findings. In scenarios where sampling coverage is too low to conduct cluster analyses, diversification rate may present an appealing alternative, although more work is needed to assess its sensitivity to sampling coverage as well. Validating these methods in alternative geographies with varying subtype distributions, sampling coverage, and transmission routes will be important to expand the generalizability of the methodology and results. The applications of lineage-level diversification rate for the study of rapidly involving pathogens, such as HIV, have only begun to be explored. 83 References 1. World Health Organization. Global Health Estimates 2016: Deaths by Cause, Age, Sex, by Country and by Region, 2000-2016. 1–8 (2018). Available at: https://www.who.int/news-room/fact-sheets/detail/the-top-10-causes-of-death. (Accessed: 23rd January 2019) 2. World Health Organization. World health statistics 2018: Monitoring health for the SDGs, sustainable development goals. 1–100 (2018). 3. World Health Organization. Ten threats to global health in 2019. 1–17 (2019). Available at: https://www.who.int/emergencies/ten-threats-to-global-health-in-2019. (Accessed: 23rd January 2019) 4. Leitner, T., Escanilla, D., Franzen, C., Uhlen, M. & Albert, J. Accurate reconstruction of a known HIV-1 transmission history by phylogenetic tree analysis. Proc Natl Acad Sci USA 93, 10864–10869 (1996). 5. Leigh Brown, A. J. et al. Transmission Network Parameters Estimated From HIV Sequences for a Nationwide Epidemic. J Infect Dis 204, 1463–1469 (2011). 6. Poon, A., Gustafson, R., Daly, P., Zerr, L. & Demlow, S. E. Near real-time monitoring of HIV transmission hotspots from routine HIV genotyping: an implementation case study. The Lancet HIV 3, e231–e238 (2016). 7. Little, S. J. et al. Using HIV Networks to Inform Real Time Prevention Interventions. PLoS ONE 9, e98443–8 (2014). 8. Centers for Disease Control CDC. Kaposi's sarcoma and Pneumocystis pneumonia among homosexual men--New York City and California. MMWR Morb. Mortal. Wkly. Rep. 30, 305–308 (1981). 9. Haverkos, H. & Curran, J. The current outbreak of Kaposi's sarcoma and opportunistic infections. CA 32, 330–339 (1982). 10. Marmor, M., Laubenstein, L., & William, D. C. Risk factors for Kaposi's sarcoma in homosexual men. The Lancet 319, 1083–1087 (1982). 11. Boyko, W. J. et al. The Vancouver Lymphadenopathy-AIDS Study: 5. Antecedent behavioural, clinical and laboratory findings in patients with AIDS and HIV-seropositive controls. CMAJ 135, 881–887 (1986). 84 12. Barré-Sinoussi, F. et al. Isolation of a T-lymphotropic retrovirus from a patient at risk for acquired immune deficiency syndrome (AIDS). Science 220, 868–871 (1983). 13. International Committee on the Taxonomy of Viruses. What to call the AIDS virus? Nature 321, 10–10 (1986). 14. Wertheim, J. O. & Worobey, M. Dating the age of the SIV lineages that gave rise to HIV-1 and HIV-2. PLoS Comput Biol 5, e1000377 (2009). 15. Silvestri, G. et al. Nonpathogenic SIV infection of sooty mangabeys is characterized by limited bystander immunopathology despite chronic high-level viremia. Immunity 18, 441–452 (2003). 16. Gao, F. et al. Origin of HIV-1 in the chimpanzee Pan troglodytes troglodytes. Nature 397, 436–441 (1999). 17. Hirsch, V. M., Olmsted, R. A., Murphey-Corb, M., Purcell, R. H. & Johnson, P. R. An African primate lentivirus (SIVsm) closely related to HIV-2. Nature 339, 389-392 (1989). 18. Visseaux, B., Damond, F., Matheron, S., Descamps, D. & Charpentier, C. HIV-2 molecular epidemiology. Infect. Genet. Evol. 46, 233–240 (2016). 19. Worobey, M. et al. Direct evidence of extensive diversity of HIV-1 in Kinshasa by 1960. Nature 455, 661–664 (2008). 20. Keele, B. F. et al. Chimpanzee reservoirs of pandemic and nonpandemic HIV-1. Science 313, 523–526 (2006). 21. Hemelaar, J., Gouws, E., Ghys, P. D. & Osmanov, S. Global and regional distribution of HIV-1 genetic subtypes and recombinants in 2004. AIDS 20, W13–W23 (2006). 22. Poon, A. F. Y. et al. The impact of clinical, demographic and risk factors on rates of HIV transmission: a population-based phylogenetic analysis in British Columbia, Canada. J. Infect. Dis. 211, 926–935 (2015). 23. UNAIDS. Fact sheet - Latest global and regional statistics on the status of the AIDS epidemic. 1–6 (2018). 24. UNAIDS. UNAIDS Fact Sheet November 2016. 1–8 (2016). 25. Mathers, B. M. et al. Global epidemiology of injecting drug use and HIV among people who inject drugs: a systematic review. The Lancet 372, 1733–1745 (2008). 26. Valera, P., Chang, Y. & Lian, Z. HIV risk inside U.S. prisons: a systematic review of risk reduction interventions conducted in U.S. prisons. AIDS Care 29, 943–952 (2016). 85 27. Schechter, M. T. et al. The Vancouver Lymphadenopathy - Aids Study: 6. HIV Seroconversion in a Cohort of Homosexual Men. CMAJ 135, 1355–1360 (1986). 28. Tyndall, M. W. et al. Intensive injection cocaine use as the primary risk factor in the Vancouver HIV-1 epidemic. AIDS 17, 887–893 (2003). 29. Hyshka, E., Strathdee, S., Wood, E. & Kerr, T. Needle exchange and the HIV epidemic in Vancouver: Lessons learned from 15 years of research. International Journal of Drug Policy 23, 261–270 (2012). 30. Montaner, J. S. G. et al. Expansion of HAART coverage is associated with sustained decreases in HIV/AIDS morbidity, mortality and HIV transmission: the ‘HIV Treatment as Prevention’ experience in a Canadian setting. PLoS ONE 9, e87872 (2014). 31. BC Centre for Disease Control. HIV in British Columbia: Annual Surveillance Report 2016. Retrieved from http://www.bccdc.ca/search?k=hiv%20annual%20report. 1–46 (2018). 32. Pinkerton, S. D. How many HIV infections are prevented by Vancouver Canada's supervised injection facility? International Journal of Drug Policy 22, 179–183 (2011). 33. Tyndall, M. W. et al. HIV seroprevalence among participants at a Supervised Injection Facility in Vancouver, Canada: implications for prevention, care and treatment. Harm Reduct J 3, 36 (2006). 34. Landovitz, R. J. et al. Epidemiology, Sexual Risk Behavior, and HIV Prevention Practices of Men who Have Sex with Men Using GRINDR in Los Angeles, California. J Urban Health 90, 729–739 (2012). 35. Lima, V. D. et al. Progress towards the United Nations 90-90-90 and 95-95-95 targets: the experience in British Columbia, Canada. J Int AIDS Soc 20, e25011 (2017). 36. BC Centre for Excellence in HIV/AIDS. About Us. 1–2 (2019). Available at: http://cfenet.ubc.ca/about-us. (Accessed: 1st February 2019) 37. Mitsuya, H. et al. 3‘-Azido-3’-deoxythymidine (BW A509U): an antiviral agent that inhibits the infectivity and cytopathic effect of human T-lymphotropic virus type III/lymphadenopathy-associated virus in vitro. Proc Natl Acad Sci USA 82, 7096–7100 (1985). 38. Fischl, M. A. et al. A randomized controlled trial of a reduced daily dose of zidovudine in patients with the acquired immunodeficiency syndrome. The AIDS Clinical Trials Group. N Engl J Med 323, 1009–1014 (1990). 86 39. Iyidogan, P. & Anderson, K. S. Current perspectives on HIV-1 antiretroviral drug resistance. Viruses 6, 4095–4139 (2014). 40. Hammer, S. M. et al. A trial comparing nucleoside monotherapy with combination therapy in HIV-infected adults with CD4 cell counts from 200 to 500 per cubic millimeter. Mass Medical Soc 335, 1081–1090 (1996). 41. Montaner, J. S. G., Hogg, R. S. & O'Shaughnessy, M. V. Emerging international consensus for use of antiretroviral therapy. The Lancet 349, 1042–1 (1997). 42. Hogg, R. S. et al. Antiviral effect of double and triple drug combinations amongst HIV-infected adults: lessons from the implementation of viral load-driven antiretroviral therapy. AIDS 12, 279–284 (1998). 43. Vernazza, P. L. et al. Effect of antiviral treatment on the shedding of HIV-1 in semen. AIDS 11, 1249–1254 (1997). 44. Cu-Uvin, S. et al. Effect of highly active antiretroviral therapy on cervicovaginal HIV-1 RNA. AIDS 14, 415–421 (2000). 45. Avettand-Fenoel, V. et al. HIV-DNA in rectal cells is well correlated with HIV-DNA in blood in different groups of patients, including long-term non-progressors. AIDS 22, 1880–1882 (2008). 46. De Cock, K. M. et al. Prevention of mother-to-child HIV transmission in resource-poor countries: translating research into policy and practice. JAMA 283, 1175–1182 (2000). 47. Quinn, T. C. et al. Viral load and heterosexual transmission of human immunodeficiency virus type 1. Rakai Project Study Group. N Engl J Med 342, 921–929 (2000). 48. Wood, E., Milloy, M. J. & Montaner, J. S. G. HIV treatment as prevention among injection drug users. Current Opinion in HIV and AIDS 7, 151–156 (2012). 49. Montaner, J. S. G. et al. The case for expanding access to highly active antiretroviral therapy to curb the growth of the HIV epidemic. The Lancet 368, 531–536 (2006). 50. Eisinger, R. W., Dieffenbach, C. W. & Fauci, A. S. HIV Viral Load and Transmissibility of HIV Infection. JAMA 321, 451–2 (2019). 51. Montaner, J. S. G. et al. Association of highly active antiretroviral therapy coverage, population viral load, and yearly new HIV diagnoses in British Columbia, Canada: a population-based study. The Lancet 376, 532–539 (2010). 87 52. World Health Organization. Guidance on preexposure oral prophylaxis (PrEP) for serodiscordant couples, men and transgender women who have sex with men at high risk of HIV. 1–21 (2012). 53. Fonner, V. A. et al. Effectiveness and safety of oral HIV preexposure prophylaxis for all populations. AIDS 30, 1973–1983 (2016). 54. Meloni, S. T. et al. Implication of First-Line Antiretroviral Therapy Choice on Second-Line Options. Open Forum Infectious Diseases 4, 1127–9 (2017). 55. Lima, V. D. et al. Differential impact of adherence on long-term treatment response among naive HIV-infected individuals. AIDS 22, 2371–2380 (2008). 56. Lima, V. D. et al. The Combined Effect of Modern Highly Active Antiretroviral Therapy Regimens and Adherence on Mortality Over Time. J. Acquir. Immune Defic. Syndr. 50, 529–536 (2009). 57. Hiener, B. et al. Identification of Genetically Intact HIV-1 Proviruses in Specific CD4+ T Cells from Effectively Treated Participants. CellReports 21, 813–822 (2017). 58. Jones, B. R. et al. Phylogenetic approach to recover integration dates of latent HIV sequences within-host. Proc Natl Acad Sci USA 115, E8958–E8967 (2018). 59. Harrigan, P. R. et al. Predictors of HIV drug-resistance mutations in a large antiretroviral-naive cohort initiating triple antiretroviral therapy. J INFECT DIS 191, 339–347 (2005). 60. Rhee, S.-Y. et al. HIV-1 Protease, Reverse Transcriptase, and Integrase Variation. Journal of Virology 90, 6058–6070 (2016). 61. Brumme, Z. L. et al. HLA-Associated Immune Escape Pathways in HIV-1 Subtype B Gag, Pol and Nef Proteins. PLoS ONE 4, e6687–13 (2009). 62. Brumme, Z. L. et al. Extensive host immune adaptation in a concentrated North American HIV epidemic. AIDS 32,1–12 (2018). 63. Kim, Y., Anderson, J. L. & Lewin, S. R. Getting the “Kill” into “Shock and Kill”: Strategies to Eliminate Latent HIV. Cell Host & Microbe 23, 14–26 (2018). 64. Vrancken, B. et al. The Genealogical Population Dynamics of HIV-1 in a Large Transmission Chain: Bridging within and among Host Evolutionary Rates. PLoS Comput Biol 10, e1003505–16 (2014). 65. Goodwin, S., McPherson, J. D. & McCombie, W. R. Coming of age: ten years of next- generation sequencing technologies. Nature Publishing Group 17, 333–351 (2016). 88 66. Goldstein, S., Beka, L., Graf, J. & Klassen, J. L. Evaluation of strategies for the assembly of diverse bacterial genomes using MinION long-read sequencing. BMC Genomics 20, 1–17 (2019). 67. Jain, M. et al. Nanopore sequencing and assembly of a human genome with ultra-long reads. Nat Biotechnol 36, 338–345 (2018). 68. Kafetzopoulou, L. E. et al. Metagenomic sequencing at the epicenter of the Nigeria 2018 Lassa fever outbreak. Science 363, 74–77 (2019). 69. Sanderson, M. J., McMahon, M. M. & Steel, M. Terraces in phylogenetic tree space. Science 333, 448–450 (2011). 70. Drummond, A. J. & Rambaut, A. BEAST: Bayesian evolutionary analysis by sampling trees. BMC Evol Biol 7, 214–8 (2007). 71. Möller, S., Plessis, du, L., the, T. S. P. O.2018. Impact of the tree prior on estimating clock rates during epidemic outbreaks. National Acad Science doi:10.1073/pnas.1713314115/-/DCSupplemental 72. Barido-Sottani, J. et al. Taming the BEAST—A Community Teaching Material Resource for BEAST 2. Syst Biol 1–5 (2017). 73. Stamatakis, A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics 30, 1312–1313 (2014). 74. Price, M. N., Dehal, P. S. & Arkin, A. P. FastTree 2--approximately maximum-likelihood trees for large alignments. PLoS ONE 5, e9490 (2010). 75. Wertheim, J. O. et al. Growth of HIV-1 Molecular Transmission Clusters in New York City. J INFECT DIS 5, e1000590–11 (2018). 76. Villandre, L. et al. Assessment of Overlap of Phylogenetic Transmission Clusters and Communities in Simple Sexual Contact Networks: Applications to HIV-1. PLoS ONE 11, e0148459–18 (2016). 77. Barido-Sottani, J., Vaughan, T. G. & Stadler, T. Detection of HIV transmission clusters from phylogenetic trees using a multi-state birth–death model. J R Soc Interface 15, 1–12 (2018). 78. Wertheim, J. O. et al. Social and Genetic Networks of HIV-1 Transmission in New York City. PLoS Pathog 13, e1006000–19 (2017). 89 79. Hué, S., Clewley, J. P., Cane, P. A. & Pillay, D. HIV-1 pol gene variation is sufficient for reconstruction of transmissions in the era of antiretroviral therapy. AIDS 18, 719–728 (2004). 80. Ragonnet-Cronin, M. et al. Automated analysis of phylogenetic clusters. BMC Bioinformatics 14, 1–1 (2013). 81. Poon, A. F. Y. Impacts and shortcomings of genetic clustering methods for infectious disease outbreaks. Virus Evolution 2, vew031–9 (2016). 82. McCloskey, R. M. & Poon, A. F. Y. A model-based clustering method to detect infectious disease transmission outbreaks from sequence variation. PLoS Comput Biol 13, e1005868–17 (2017). 83. Stadler, T. & Bonhoeffer, S. Uncovering epidemiological dynamics in heterogeneous host populations using phylogenetic methods. Philos. Trans. R. Soc. Lond., B, Biol. Sci. 368, 20120198 (2013). 84. Hall, M. D., Woolhouse, M. E. J. & Rambaut, A. The effects of sampling strategy on the quality of reconstruction of viral population dynamics using Bayesian skyline family coalescent methods: A simulation study. Virus Evolution 2, 1–14 (2016). 85. Coltart, C. E. M. et al. Ethical considerations in global HIV phylogenetic research. The Lancet HIV 1–11 (2018). doi:10.1016/S2352-3018(18)30134-6 86. Romero-Severson, E., Skar, H., Bulla, I., Albert, J. & Leitner, T. Timing and order of transmission events is not directly reflected in a pathogen phylogeny. Molecular Biology and Evolution 31, 2472–2482 (2014). 87. Romero-Severson, E. O., Bulla, I. & Leitner, T. Phylogenetically resolving epidemiologic linkage. Proc Natl Acad Sci USA 113, 2690–2695 (2016). 88. Bernard, E. J., Azad, Y., Vandamme, A. M., Weait, M. & Geretti, A. M. HIV forensics: pitfalls and acceptable standards in the use of phylogenetic analysis as evidence in criminal investigations of HIV transmission. HIV Med 8, 382–387 (2007). 89. Gilbert, M., Swenson, L., Unger, D., Scheim, A. & Grace, D. Need for robust and inclusive public health ethics review of the monitoring of HIV phylogenetic clusters for HIV prevention. The Lancet HIV 3, e461 (2016). 90. Canadian HIV/AIDS Legal Network. Criminal Law & HIV Non-Disclosure in Canada. 1–18 (2014). 90 91. Ou, C. Y. et al. Molecular Epidemiology of HIV Transmission in a Dental Practice. Science 256, 1165–1171 (1992). 92. Scaduto, D. I. et al. Source identification in two criminal cases using phylogenetic analysis of HIV-1 DNA sequences. Proced National Acad Sciences 50, 21242-21247 (2010). 93. Council for International Organizations of Medical Sciences World Health Organization. International Ethical Guidelines for Health-related Research Involving Humans. 1–119 (2017). Available at: https://cioms.ch/wp-content/uploads/2017/01/WEB-CIOMS-EthicalGuidelines.pdf. (Accessed: 2nd March 2019) 94. Jetz, W., Thomas, G. H., Joy, J. B., Hartmann, K. & Mooers, A. O. The global diversity of birds in space and time. Nature 491, 444–448 (2012). 95. Joy, J. B. et al. Phylogenetically estimated HIV diversification rates reveal prevention of HIV-1 by antiretroviral therapy. Journal of the International AIDS Society 16, 1–1 (2015). 96. Faith, D. Conservation evaluation and phylogenetic diversity. The Lancet 61, 1–10 (1992). 97. Redding, D. W. et al. Evolutionarily distinctive species often capture more phylogenetic diversity than expected. Journal of Theoretical Biology 251, 606–615 (2008). 98. Magallón, S. & Sanderson, M. Absolute diversification rates in angiosperm clades. BioOne 55, 1762 (2001). 99. Vane-Wright, R. I., Humphries, C. J. & Williams, P. H. What to protect?—Systematics and the agony of choice. Biological Conservation 55, 235–254 (1991). 100. Pavoine, S., Ollier, S. & Dufour, A. B. Is the originality of a species measurable? Ecology Letters 8, 579–586 (2005). 101. May, R. M. Taxonomy as destiny. Nature 347, 129–130 (1990). 102. Altschul, S. F. & Lipman, D. J. Equal Animals. Nature 348, 493–494 (1990). 103. Isaac, N. J. B., Turvey, S. T., Collen, B., Waterman, C. & Baillie, J. E. M. Mammals on the EDGE: Conservation Priorities Based on Threat and Phylogeny. PLoS ONE 2, e296–7 (2007). 104. Redding, D. W. & Mooers, A. O. Incorporating evolutionary measures into conservation prioritization. Conserv. Biol. 20, 1670–1678 (2006). 91 105. Cadotte, M. W. et al. Phylogenetic diversity metrics for ecological communities: integrating species richness, abundance and evolutionary history. Ecology Letters 13, 96–105 (2010). 106. Rao, C. R. Diversity and Dissimilarity Coefficients - a Unified Approach. Theoretical Population Biology 21, 24–43 (1982). 107. Pauplin, Y. Direct Calculation of a Tree Length Using a Distance Matrix. J Mol Evol 51, 41–47 (2000). 108. Gudde, R. M., Joy, J. B. & Mooers, A. O. Imperilled phylogenetic endemism of Malagasy lemuriformes. Diversity and Distributions 19, 664–675 (2013). 109. Guerin, G. R., Ruokolainen, L. & Lowe, A. J. A georeferenced implementation of weighted endemism. Methods Ecol Evol 6, 845–852 (2015). 110. Guerin, G. R. & Lowe, A. J. Mapping phylogenetic endemism in R using georeferenced branch extents. SoftwareX 3-4, 22–26 (2015). 111. Ragonnet-Cronin, M., of, S. L. T. J.2015. Transmission of non-B HIV subtypes in the United Kingdom is increasingly driven by large non-heterosexual transmission clusters. academic.oup.com doi:10.1093/infdis/jiv758/-/DC1 112. Ratmann, O. et al. Phylogenetic Tools for Generalized HIV-1 Epidemics: Findings from the PANGEA-HIV Methods Comparison. Molecular Biology and Evolution 34, 185–203 (2017). 113. Holmes, E. C., Nee, S., Rambaut, A., Garnett, G. P. & Harvey, P. H. Revealing the history of infectious disease epidemics through phylogenetic trees. Philosophical Transactions of the Royal Society B: Biological Sciences 349, 33–40 (1995). 114. Kosakovsky Pond, S. L. et al. An Evolutionary Model-Based Algorithm for Accurate Phylogenetic Breakpoint Mapping and Subtype Prediction in HIV-1. PLoS Comput Biol 5, e1000581–21 (2009). 115. Katoh, K. & Toh, H. Recent developments in the MAFFT multiple sequence alignment program. Briefings in Bioinformatics 9, 286–298 (2008). 116. Larsson, A. AliView: a fast and lightweight alignment viewer and editor for large datasets. Bioinformatics 30, 3276–3278 (2014). 117. Rhee, S.-Y. et al. Geographic and Temporal Trends in the Molecular Epidemiology and Genetic Mechanisms of Transmitted HIV-1 Drug Resistance: An Individual-Patient- and Sequence-Level Meta-Analysis. PLoS Med 12, e1001810–29 (2015). 92 118. Paradis, E., Claude, J. & Strimmer, K. APE: Analyses of Phylogenetics and Evolution in R language. Bioinformatics 20, 289–290 (2004). 119. Rambaut, A., Lam, T. T., Max Carvalho, L. & Pybus, O. G. Exploring the temporal structure of heterochronous sequences using TempEst (formerly Path-O-Gen). Virus Evolution 2, vew007–7 (2016). 120. Statistics Canada. Dictionary, Census of Population, 2016 - Census tract (CT). 1–2 (2016). Available at: https://www12.statcan.gc.ca/census-recensement/2016/ref/dict/geo013-eng.cfm. (Accessed: 3rd January 2019) 121. Hu, M.-C., Pavlicova, M. & Nunes, E. V. Zero-Inflated and Hurdle Models of Count Data with Extra Zeros: Examples from an HIV-Risk Reduction Intervention Trial. The American Journal of Drug and Alcohol Abuse 37, 367–375 (2011). 122. Wood, E. et al. Longitudinal community plasma HIV-1 RNA concentrations and incidence of HIV-1 among injecting drug users: prospective cohort study. BMJ 338, b1649 (2009). 123. Xia, Q., Wiewel, E. W. & Torian, L. V. Revisiting the methodology of measuring HIV community viral load. J. Acquir. Immune Defic. Syndr. 63, e82–4 (2013). 124. Dark, S. J. & Bram, D. The modifiable areal unit problem (MAUP) in physical geography. Progress in Physical Geography 31, 471–479 (2016). 125. Wickham, H., Francois, R., Henry, L. & Müller, K. dplyr: A Grammar of Data Manipulation. R package version 0.7.4. https://CRAN.R-project.org/package=dplyr (2016). 126. Wickham, H. & Henry, L. tidyr: Easily Tidy Data with 'spread()' and ‘gather()’ Functions. R package version 0.8.0. https://CRAN.R-project.org/package=tidyr. 1–22 (2018). 127. Bivand, R., Keitt, T. & Rowlingson, B. rgdal: Bindings for the ‘Geospatial’ Data Abstraction Library. R package version 1.2-16. https://CRAN.R-project.org/package=rgdal. (2017). 128. Bivand, R. S. & Lewin-Koh, N. maptools: Tools for Handling Spatial Objects. R package version 0.9-4. https://CRAN.R-project.org/package=maptools. 1–83 (2018). 129. Hijmans, R. J. raster: Geographic Data Analysis and Modeling. R package version 2.6-7. https://CRAN.R-project.org/package=raster. (2017). 93 130. Baquero, O. S. ggsn: North Symbols and Scale Bars for Maps Created with 'ggplot2' or ‘ggmap’. R package version 0.4.0. https://CRAN.R-project.org/package=ggsn. 1–8 (2017). 131. Brunsdon, C. & Chen, H. GISTools: Some further GIS capabilities for R. R package version 0.7-4. https://CRAN.R-project.org/package=GISTools. 1–25 (2014). 132. Bivand, R. S., Pebesma, E. & Gomez-Rubio, V. Classes and methods for spatial data in R. R News 5 (2), https://cran.r-project.org/doc/Rnews/. 1–118 (2005). 133. Peng, R. D., Murdoch, D., Rowlingson, B. & Murta, A. gpclib: General Polygon Clipping Library for R. R package version 1.5-5. https://CRAN.R-project.org/package=gpclib. 1–9 (2013). 134. Bivand, R. & Rundel, C. rgeos: Interface to Geometry Engine - Open Source (‘GEOS’). R package version 0.4-1. https://CRAN.R-project.org/package=rgeos. 1–78 (2018). 135. Wickham, H. ggplot2: Elegant Graphics for Data Analysis. 1–222 (2016). 136. BC Centre for Disease Control. HIV in British Columbia: Annual Surveillance Report 2015. Retrieved from http://www.bccdc.ca/search?k=hiv%20annual%20report. 1–45 (2017). 137. To, T. H., Jung, M., Lycett, S. & Gascuel, O. Fast dating using least-squares criteria and algorithms. Syst Biol 65, 82-97 (2016). 138. Lythgoe, K. A., Gardner, A., Pybus, O. G. & Grove, J. Short-Sighted Virus Evolution and a Germline Hypothesis for Chronic Viral Infections. Trends in Microbiology 25, 336–348 (2017). 139. Hallett, T. B., White, P. J. & Garnett, G. P. Appropriate evaluation of HIV prevention interventions: from experiment to full-scale implementation. Sexually Transmitted Infections 83, i55–i60 (2007). 140. Smith, M. K., Powers, K. A., Muessig, K. E., Miller, W. C. & Cohen, M. S. HIV Treatment as Prevention: The Utility and Limitations of Ecological Observation. PLoS Med 9, e1001260–9 (2012). 94 Appendices Appendix A Supplementary material for Chapter 4: Viral diversification rate identifies HIV transmission hot spots A.1 Supplementary methods Sequences The clinical and research staff at the BC-CfE conducted all laboratory work including RNA extraction, PCR, and sequencing. The majority of partial pol sequences were 1,497 bases (b) in length (n=30,445), encoding protease and the first 400 codons of reverse transcriptase. Another group of sequences (n=5,118) spanned 1,017 b, which covered all of protease and the first 240 codons of reverse transcriptase. Finally, a group of 709 sequences were assembled from three partial sequences of 276 b, 297 b, and 363 b, for a total of a 936 b concatenated sequence covering protease and reverse transcriptase codons 24-236. Select samples additionally had the entire integrase gene sequenced (n=3,867, representing 2,024 unique PLHIV) covering 864 b. Integrase sequences were analyzed in the same manner as the partial pol sequences as a validation of inferred trees and resulting diversification rates (Table A-1). Clinical, epidemiological, and demographic variables In addition to sequences, available data for this analysis included sample collection date; treatment initiation date; treatment regimen; date of first viral load; physician estimated date of seroconversion; plasma viral load (HIV RNA copies/mL); HIV subtype classifications (generated using the SCUEAL algorithm114); ethnicity; birth year; sex at birth; sex; self-reported risk factors (injection drug use, men who have sex with men, heterosexual contact, any receipt of blood product 95 or exposure to blood risk, other risk exposure); having ever tested positive for hepatitis C infection; having ever had acquired immune deficiency syndrome (AIDS); if applicable date of mortality and cause of death; local health authority, health authority, postal code, census tract (CT), and census metropolitan area (CMA) of patient residence; and forward sortation area (first three digits of postal code) of physician requesting test. The viral load assay detection limit changed several times prior to 2008, so for consistency, only data from 2008 onwards was analyzed in the final model. The detection limit over that period was 40–10,000,000 copies viral RNA/mL plasma using a Roche COBAS HIV-1 Ampliprep Taqman assay51. To suit the longitudinal nature of the analysis, patient data was carried over year to year until new data was available (ie. more recent viral load), unless participants were removed or passed away. HIV sequences and patient information were stored in secure, access-restricted facilities at the BC-CfE in a secure, password protected Oracle database. Patient data were de-identified and doubly anonymized with randomly generated 6-character identifiers. No document was created linking identifiers to the patient. Diversification rate robustness In contrast to clade-level measures of diversification98, lineage-level diversification rate is an integrated splitting rate that is unique to each tip in the tree and considers the entire path from root to tip94. It is the reciprocal of the equal splits metric of evolutionary distinctiveness, which divides the phylogenetic distance of each edge by its daughter edges until the entire tree has been accounted for104. Since equal splits is an indication of tips’ evolutionary uniqueness, its reciprocal, diversification rate, is an estimate of the extent to which tips are evolutionarily shared. A virus in a rapidly diversifying clade is typified by short branch lengths shared among many tree neighbours (Figure A-4), resulting in a low equal splits value and a high diversification rate94. Since the viral phylogenetic tree is an estimate of a transmission tree, the lineage-level diversification rate can be 96 interpreted as an estimate of the between-host transmission rate for epidemiological purposes. If an individual’s diversification rate increases over the year, either they were directly involved in a recent transmission event or one or more individuals in close phylogenetic proximity to them were involved in a recent transmission event. A binary parameter to detect changes in individuals’ viral diversification rate greater than 1 (subs/site)-1 was calculated to bring greater focus to recent and proximate transmission events. The diversification rates obtained through an approximate maximum likelihood with root-to-tip regression framework were compared to those obtained from time-scaled trees (implemented in Least Squares Dating), maximum likelihood trees, outgroup-rooted approximate maximum likelihood trees, and integrase approximate maximum likelihood trees (Table A-1). Aggregated variables Patient-level variables were aggregated according to the merged census tract of patients’ residence using R packages dplyr125 and tidyr126. For each continuous quantitative variable - including viral load, diversification rate, change in diversification rate, and age - the sum, mean, median, maximum, mean of the top 3, mean of the top 5 were calculated. For discrete count variables – including prevalent cases, new cases, and binary change in diversification rate – the sum was calculated. For all sums, the population normalized sums (divided by the total population of the merged census tract) and area-normalized sum (divided by the area in km2 of the merged census tract) were also calculated. The census tract population from the nearest census (census years within study period 2006, 2011, 2016) was used. For study years with equidistant census years, the average population between the two populations was taken. Qualitative variables – including sex, sex at birth, risk factors, and ethnicities – were aggregated as sums and then further evaluated 97 as proportions of the total population of PLHIV in a given merged census tract. The total number of PLHIV in a merged census tract reflects the number of people connected to care with sequences available through the Drug Treatment Program living in that census tract in that year. Individuals who have ever had a detectable viral load measured in BC were carried forward to subsequent time intervals until they were removed from the cohort due to death or migration out of province. There was extensive collinearity among some aggregated predictor variables, which was expected given how many aggregate variables were sums across the prevalent population. Collinearity was taken into account when evaluating variables for inclusion in the final model (Figure A-7). The outcome, the total number of new HIV cases within each merged census tract in the subsequent year, was estimated based on the patients’ date of first detectable viral load. No further imputation of the infection date or the unsampled population was conducted in part because the excessive parameterization required for most imputation methods is not robust for the small populations in merged census tracts. The population-normalized incidence and area-normalized incidence were also evaluated as potentially informative outcomes, but the total new cases outcome was deemed more informative towards the burden on the health care system. Our overall provincial estimates of cumulative annual incidence were compared to the BC Centre for Disease Control (BCCDC) estimates using reported new diagnoses, as well as BCCDC estimates of true incidence (Figure A-1). Spatial analyses An ESRI shapefile (A.3 Supplementary data and shapefiles) delineating census tracts in Canada was obtained from the 2016 census on the Statistics Canada website as a Lambert conformal conic projected of North American Datum of 1983 (NAD83), which was then converted to un-projected 98 latitude/longitude coordinates for visualization (Figure A-5). The population and area in square kilometers of each census tract for all recent censuses (1996, 2001, 2006, 2011, and 2016) were also obtained from Statistics Canada (A.3 Supplementary data and shapefiles). The shapefiles were read into R using the rgdal127 package. Other R packages that were employed for spatial analyses included maptools128, raster129, ggsn130, GIStools131, sp132, gpclib133, and rgeos134. The census tract polygons were united to suit the study in two ways. Firstly, to account for changes in census tract borders over the course of the study due to population changes, the decimals of the census tract were truncated off and all census tracts sharing a common ID were united. Secondly, in order to maintain patient confidentiality, census tracts were iteratively united based on proximity and sharing the same census metropolitan area with the goal of every polygon containing a minimum of three PLHIV in all of the time intervals (Figure A-5). The final polygons are referred to as merged census tracts or census tracts for simplicity. After building the trees and calculating diversification rate, samples from patients with no recorded census tract were removed from the analysis (n=2,686). The Statistics Canada postal code conversion file (A.3 Supplementary data and shapefiles) was used to generate census tracts for patients with no recorded census tract who had a recorded postal code (n=243). Chloropleth heat maps of merged census tract polygons were plotted using a combination of dplyr125, tidyr126, ggplot2135, and RColorBrewer34. Model training and testing A simple supervised machine learning method was applied to identify a predictive model, specifically by using a holdout cross-validation method (Figure A-7). The data was segregated into a training subset (2009, 2010, 2012, 2013, 2015, and 2016) and a testing subset of interspersed years (2008, 2011, and 2014). The 2017 dataset (2017 predictors, 2018 new cases) was reserved 99 for a final blind prediction and the full dataset for 2018 was obtained after the final model was selected. Bivariate relationships between all variables were calculated as Pearson’s correlation coefficients to inform their prioritized inclusion in the model as well as to assess potential collinearity. Since ZINBs are two-part models, each variable was evaluated separately in the binomial and count component of the model for inclusion, resulting in some variables only being included in one part of the model. In order for a variable to be included, the full model had to significantly reduce the |log-likelihood| relative to its nested version, as assessed by a likelihood ratio test with a significance level of 0.05; the full model had to have a lower AIC and BIC; and the coefficient p-value must also have been <0.05. A.2 Supplementary results and discussion A.2.1 Outcome validity Our estimates of new HIV cases in BC based on the date of first detectable viral load from the DTP cohort with sequences were consistent with the BCCDC estimates over the study period (Figure A-1). For instance, in 2014, the BCCDC estimated the total number of new HIV infections was 305 (95% CI, 210-400 cases)136; our estimate was within that range at 293 estimated new cases in 2014. This difference could be partially attributable to the fact that sometimes initial sequences are generated without the new case being reported to the CDC, as would be the case for a new migrant to BC. Although the overall provincial incidence steadily decreased since 2008 down to 5.1 cases per 100,000 in 201631, this does not reflect the persistently high incidence among several communities. Notably, we observed a maximum census tract-level HIV incidence of 75.4 cases per 100,000 in 2017 within a central census tract in the West End neighbourhood of Vancouver in close proximity to St. Paul’s hospital where this study was conducted. 100 While patients’ date of first detectable viral load is not equivalent to date of infection, it is a suitable estimate for this purpose. Firstly, since it is used to estimate cumulative new cases over the course of the year, it is flexible to some uncertainty in the exact date of infection within that year. Further, current methods of estimating date of infection rely on highly parametrized models of CD4 counts, for example, that would introduce a great deal of uncertainty into the model and which were unlinked to the patient data for this study. Although it is possible that a decrease in new HIV cases could be the result of fewer people getting tested as opposed to fewer new infections, this is likely not the case in BC where the percent of PLHIV who have been diagnosed is estimated to be increasing over time35. A.2.2 Diversification rate robustness and characterization The calculation of individual diversification rates was robust to the tree inference software, the rooting method, and the scale of the tree (Table A-1). We validated that our derivation of diversification rates was not affected by the tree rooting method by comparing outgroup rooting to root-to-tip regression rooting across the full set of trees. The trees were virtually identical as were the diversification rates, resulting in a mean correlation of 1.00. Further, we validated that diversification rates calculated from time-scaled trees generated using least squares dating137 correlated well with the standard subs/site trees, with a mean Pearson correlation of 0.582 and a mean Spearman correlation of 0.817. Lastly, we found that the standard partial pol trees led to similar diversification rates as those calculated from the integrase gene (mean Pearson correlation, 0.434; mean Spearman correlation, 0.559). The correlation improves over time as the number of patients with both partial pol and integrase increases. Our estimates of patient-level and 101 geographically-aggregated diversification rates were robust to different tree-building methods, consistent with Joy et al.94, however some variation in tree topology is inherent when analyzing large distributions of phylogenies with many tips. It is worth noting that we only used the oldest patient sample available at a given time interval to infer trees and to derive diversification rates. The inclusion of more than one sample per patient per tree would artificially inflate the diversification rates for patients with the most samples. We justified using the oldest available sequence because founder theory suggests that early HIV variants within a host are better adapted to transmission than later variants, which tend to be adapted for immune evasion138. Because each year’s tree includes patients’ oldest sequence from that year or previously, this resulted in an increase in the size of the tree over time (Figure A-3). To detect whether there were any population-level changes in diversification rate over time, we evaluated the distribution of patients’ mean diversification rates across bootstraps and their associated coefficients of variation (standard deviation/mean*100). The distribution of each patient’s mean diversification rate across the 100 bootstrap trees for each year is approximately log normal with a right skew, characteristic of a gamma distribution (Figure A-4). The mean diversification rate across patients in each year from 2008 to 2017 varied from 98.1 (subs/site)-1 to 120.0 (subs/site)-1 while the median ranged from 62.6 (subs/site)-1 to 68.0 (subs/site)-1. A Kruskal-Wallis test revealed that the distribution of diversification rates in the entire population did not significantly change across the study period (p = 0.385). There was a slight increase in both mean and median diversification rate over time, which is attributable to sampling from a larger tree with more extensive branching. 102 A.2.3 Zero-inflated negative binomial model interpretation The relationships between new HIV cases in the subsequent year and significant predictors were summarized in Table A-3. The binomial part of the model reflects the probability that a polygon has more than zero new HIV cases and as such, the exponentiated coefficients are interpreted as adjusted odds ratio. For every additional individual added to the prevalent PLHIV population in an area, the adjusted odds of having greater than zero new HIV cases in the subsequent year was 1.23 (95% CI, 1.14-1.35) times lower. For every 1-unit increase in the mean of the top 5 ln(diversification rates) of a census tract, the adjusted odds of having greater than zero new HIV cases in the subsequent year was 3.10 (95% CI, 1.61-5.97) times higher. The negative binomial component of the model reflects the predicted number of new HIV cases, given that there are more than zero new cases, and the exponentiated coefficients are interpreted as adjusted relative risks (incident rate ratios). For every additional new HIV case in a census tract in a given year, the adjusted risk of new HIV cases in the subsequent year was 1.07 (95% CI, 1.03-1.11) times greater. For every 1-unit increase in the mean of the top 5 ln(diversification rate) of a census tract, the adjusted risk of new HIV cases in the subsequent year was 1.38 (95% CI, 1.28-1.49) times higher. For every additional PLHIV with an annual change in diversification rate ³ 1 (subs/site)-1 in a census tract, the adjusted risk of new HIV cases in the subsequent year was 1.10 (95% CI, 1.07-1.13) times higher. A.2.4 Spatial analyses and interpretations In 2016, there were 766 census tracts in BC, after truncating the decimal off of these census tracts (excluding those in the downtown core to maintain their resolution), there were 392 truncated 103 census tracts, of which 114 had only 1-2 PLHIV at some point in the study period. After merging census tracts in BC such that no polygon had fewer than 3 people living with HIV at any point in the study period, there were a total of 265 merged census tracts, 250 of which contained PLHIV at some point during the study period. Of the 2,500 estimates of number of new HIV cases within merged census tracts between 2008 and 2017, 1,468 of the observations were zero (Figure A-2), justifying the use of a zero-inflated model. When interpreting the geographic variation of any variable, it is important to consider that the way the data were aggregated can seriously impact the conclusions. This is commonly referred to as the modifiable aerial unit problem124. To ensure that the conclusions we drew were robust to changes in the borders of the study polygons, in the future we will repeat the analysis with an alternative geography, namely local health authorities. Similar neighbourhoods were identified at the top of the list of areas ranked as high risk. Further, one must be cautious in the interpretation of ecologic studies of this nature, such that a statistical relationship at one level of aggregation should not be assumed to hold at any other levels of aggregation, including at the individual level139,140. An area assigned as having a high transmission risk does not equate to an individual residing in that area as necessarily having a high probability of infection. People are not isolated in their communities and we do not exclude the reality that PLHIV from one census tract could transmit HIV to someone outside their census tract of residence. However, our work does demonstrate that many new cases arose within the same census tracts where people resided who were associated with frequent transmission. Our dataset was limited to those connected to care within BC and therefore we cannot make any statements about the interactions of our study population with people who live outside of BC. Expanding these analyses 104 to a global scale would be a worthy pursuit, although it would prove challenging to connect disparate datasets containing sequences and geography of patient residence. Further, the inter-provincial migration of PLHIV in Canada is not well tracked. Thus, it is possible that some individuals who we deemed to be newly infected based on a first viral load measurement could have been infected outside the province. Finally, we assumed that our study population did not migrate between census tracts within BC during the study period. Although this assumption is surely not true for all individuals, our dataset was restricted to the geography of patient residence at the first viral load test. 105 Figure A-1 A comparison of annual estimates of new HIV cases between the BC Centre for Disease Control and this study. (a) The total new HIV cases in British Columbia over time estimated using the date of first detectable viral load in the Drug Treatment Program cohort (red line) as compared to the BC Centre for Disease Control number of new diagnoses (blue line) and their estimates of newly infected cases with 95% confidence intervals (black). (b) The total new HIV cases in individual merged census tracts over time, as estimated using the date of first detectable viral load. 05102008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018YearNew HIV cases in merged census tracts●●2002503003504002008 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018YearTotal new HIV cases(b) (a) 106 Figure A-2 The distribution of total annual new HIV cases within merged census tracts in BC across the study period. 1,468 observations out of n= 2,500 over ten years were zero. Figure A-3 The total number of tips in each time interval phylogenetic tree. The number of tips represents the number of unique participants in the Drug Treatment Program with a detectable HIV sequence (plus one HXB2 reference sequence) between 1996 and March 2018. 80213141640216325842982340537924242470052085678614366617122758879848346867290119332 9573963102500500075001000019961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018YearTips107 Figure A-4 The distribution of patients’ mean and % coefficient of variation (CV) for diversification rates calculated across 100 bootstrap approximate maximum likelihood trees. 0.00.20.42 4 6 8ln(mean DR) of patients across bootstrapsDensityYear20082009201020112012201320142015201620170.00.20.40.60.8−6 −4 −2 0 2 4ln(%CV diversification rate) of patients across bootstrapsDensityYear2008200920102011201220132014201520162017108 Figure A-5 Merged census tracts within British Columbia, Metro Vancouver, and downtown Vancouver. Census tracts, highlighted in blue, are small geographic areas with populations generally under 10,000 people with relatively homogenous socioeconomic characteristics that are within census agglomerations with a population of 50,000 or more in the previous census106. The original 2016 census tract borders are shown in yellow and the final merged census tract borders are shown in black. 109 Figure A-6 Predictive power of the full zero-inflated negative binomial (ZINB) model compared to the reduced ZINB model over time for testing, training, and blind prediction data. Predictive criteria included the Pearson’s correlation coefficient (for observed cases > 0 to remove the effect of zero-inflation); sensitivity (true positives/ (true positives + false negatives)); specificity (true negatives/(true negatives + false positives)); positive predictive value (true positives/(true positives + false positives)); negative predictive value (true negatives/(true negatives + false negatives)), and accuracy ((true positives + true negatives)/(true positives + true negatives + false positives + false negatives)). 110 Patient variables, 2008-2017: lineage-level diversification rate + clinical, epidemiological, demographic dataAggregated variables, 2008-2017 Build Pearson correlation matrixVariables ordered by correlation with outcome• Colinear with existing variables in model (r>0.5)? • Other coefficients standard error up >10%?Null model• Xi coefficient p-value <0.05?• Lower AIC, BIC, log-likelihood vs. reduced (best) model?• Likelihood ratio test vs. reduced (best) model, p<0.05?Aggregate by patients’ merged census tract of residenceAdd variable XiTraining data: 2009, 2010, 2012, 2013, 2015, and 2016Testing data: 2008, 2011, and 2014New best modelInclude XiCo-linearity matrix+Test predictive power: sensitivity, specificity, PPV, NPV, top 10 agreement, top 20 agreementSeparate data into training, testing, and validation datasetsNoYesYesNoExclude XiFinal modelFinal validation: 2017 predictorsi = i+1i = i+1when i=total number variables evaluatedFigure A-7 Workflow of data segregation, model training, testing, and blind prediction. 111 Table A-1 A comparison of diversification rates calculated using different tree building strategies. The Pearson’s correlation coefficients and Spearman rank correlation coefficients by year comparing the partial pol diversification rates obtained from approximate maximum likelihood trees with root-to-tip regression rooting as compared to integrase trees, outgroup rooted tree, and time scaled trees. Integrase vs Partial Pol Outgroup Root vs Root-to-Tip Time-Scaled vs Subs/Site Year Pearson Spearman n Pearson Spearman n Pearson Spearman n 2008 0.192 0.407 106 1.000 1.000 6143 0.626 0.797 6143 2009 0.228 0.399 117 1.000 1.000 6661 0.613 0.800 6661 2010 0.308 0.424 130 1.000 1.000 7122 0.612 0.805 7122 2011 0.372 0.438 147 1.000 1.000 7588 0.629 0.811 7588 2012 0.429 0.458 163 1.000 1.000 7984 0.538 0.814 7984 2013 0.482 0.465 186 1.000 1.000 8346 0.512 0.818 8346 2014 0.406 0.573 254 1.000 1.000 8672 0.526 0.821 8672 2015 0.615 0.674 378 1.000 1.000 9011 0.547 0.825 9011 2016 0.650 0.758 650 1.000 1.000 9332 0.582 0.830 9332 2017 0.545 0.776 865 1.000 1.000 9573 0.614 0.834 9573 2018 0.547 0.781 913 1.000 1.000 9631 0.609 0.834 9631 Mean 0.434 0.559 - 1.000 1.000 - 0.582 0.817 - Table A-2 Pearson’s correlations between all evaluated potential aggregated predictors and the outcome. The outcome, total new HIV cases in the polygon in the subsequent year, is referred to as Sum_NextFirstVL. Variable Pearson Corr. Coeff. p-value Adj. p-value (Bonferroni) Variable Description Sum_NextFirstVL 1.00 0.00E+00 0.00E+00 The sum of new HIV cases in the subsequent year Sum_FirstVL 0.74 0.00E+00 0.00E+00 The sum of new HIV cases in a given year Sum_BinaryChangeDR 0.72 0.00E+00 0.00E+00 The sum of PLHIV with an annual change in diversification rate >=1 Sum_InfoKnown 0.71 0.00E+00 0.00E+00 The sum of PLHIV for whom all risk info was known Sum_RiskKnown 0.71 0.00E+00 0.00E+00 The sum of PLHIV for whom some risk info was known Sum_White 0.70 0.00E+00 0.00E+00 The sum of PLHIV who identify as white ethnicity Sum_log10VL 0.69 0.00E+00 0.00E+00 The sum of all log10 (viral loads) Sum_Count 0.69 0.00E+00 0.00E+00 The sum of all PLHIV Sum_Male 0.69 0.00E+00 0.00E+00 The sum of all PLHIV who are male PerTotalPop_FirstVL 0.66 3.81E-291 5.41E-289 The sum of new HIV cases in a given year per 100,000 Sum_DR 0.65 9.06E-286 1.29E-283 The sum of diversification rates 112 Sum_Other 0.64 3.95E-278 5.61E-276 The sum of PLHIV exposed to 'other' risks PerTotalPop_BinaryChangeDR 0.64 3.55E-270 5.04E-268 The sum of PLHIV with an annual change in diversification rate >=1 per 100,000 PerArea_HCV 0.63 1.09E-263 1.55E-261 The sum of PLHIV ever co-infected with HCV per km^2 PerArea_Hetero 0.63 1.25E-263 1.77E-261 The sum of PLHIV who participate in heterosexual activities per km^2 PerArea_IDU 0.63 6.16E-258 8.75E-256 The sum of PLHIV who identify as IDU per km^2 PerArea_FirstVL 0.62 1.26E-255 1.80E-253 Thes sum of new HIV cases in a given year per km^2 Sum_HCV 0.62 1.11E-246 1.57E-244 The sum of PLHIV ever co-infected with HCV PerTotalPop_InfoKnown 0.61 1.24E-245 1.76E-243 The sum of PLHIV for whom all risk info was known per 100,000 PerTotalPop_RiskKnown 0.61 1.07E-243 1.51E-241 The sum of PLHIV for whom some risk info was known per 100,000 PerTotalPop_Other 0.61 1.23E-240 1.74E-238 The sum of PLHIV exposed to 'other' risks per 100,000 Sum_MSM 0.61 1.12E-238 1.59E-236 The sum of PLHIV who identify as MSM PerTotalPop_Count 0.60 1.18E-235 1.68E-233 The sum of all PLHIV per 100,000 PerTotalPop_White 0.60 1.87E-234 2.65E-232 The sum of PLHIV who identify as white ethnicity per 100,000 PerArea_FirstNation 0.60 4.40E-233 6.25E-231 The sum of PLHIV who identify as First Nations per km^2 PerTotalPop_log10VL 0.60 6.00E-232 8.52E-230 The sum of all log10 (viral load) per 100,000 Sum_IDU 0.60 1.03E-229 1.46E-227 The sum of PLHIV who identify as IDU PerArea_Blood 0.60 2.38E-229 3.38E-227 The sum of PLHIV exposed to blood risk per km^2 PerArea_Other 0.59 3.11E-226 4.42E-224 The sum of PLHIV exposed to 'other' risks per km^2 PerArea_BinaryChangeDR 0.59 5.84E-222 8.29E-220 The sum of PLHIV with a change in diversification rate >=1 per km^2 PerTotalPop_Male 0.59 6.62E-218 9.40E-216 The sum of all PLHIV who are male per 100,000 Sum_Hetero 0.58 2.27E-212 3.22E-210 The sum of PLHIV who participate in heterosexual activities PerTotalPop_DR 0.57 1.88E-204 2.67E-202 The sum of diversification rates per 100,000 PerArea_RiskKnown 0.57 2.07E-201 2.93E-199 The sum of PLHIV for whom some risk info was known per km^2 PerArea_InfoKnown 0.56 1.53E-198 2.17E-196 The sum of PLHIV for whom all risk info was known per km^2 Sum_Blood 0.56 4.86E-198 6.90E-196 The sum of PLHIV exposed to blood risk PerArea_Count 0.55 4.67E-191 6.63E-189 The sum of all PLHIV per km^2 PerTotalPop_Blood 0.55 4.20E-188 5.97E-186 The sum of PLHIV exposed to blood risk per 100,000 Sum_FirstNation 0.55 5.61E-186 7.96E-184 The sum of PLHIV who identify as First Nations PerArea_Black 0.54 2.70E-181 3.83E-179 The sum of PLHIV who identify as black ethnicity per km^2 PerArea_log10VL 0.54 8.93E-178 1.27E-175 The sum of all log10 (viral load) per km^2 PerArea_White 0.54 1.23E-177 1.75E-175 The sum of PLHIV who identify as white ethnicity per km^2 PerTotalPop_HCV 0.53 8.09E-171 1.15E-168 The sum of PLHIV ever co-infected with HCV per 100,000 PerArea_DR 0.52 2.18E-167 3.09E-165 The sum of diversification rates per km^2 PerTotalPop_MF 0.52 1.08E-163 1.54E-161 The sum of all PLHIV who are transgender (male to female) per 100,000 113 PerArea_Male 0.52 1.19E-162 1.69E-160 The sum of all PLHIV who are male per km^2 PerTotalPop_IDU 0.51 6.87E-158 9.76E-156 The sum of PLHIV who identify as IDU per 100,000 PerArea_Female 0.51 6.32E-156 8.97E-154 The sum of all PLHIV who are female per km^2 PerTotalPop_MSM 0.51 5.47E-155 7.77E-153 The sum of PLHIV who identify as MSM per 100,000 Sum_changeDR 0.50 2.80E-152 3.97E-150 The sum of all annual changes in diversification rate PerArea_MF 0.50 5.97E-150 8.48E-148 The sum of all PLHIV who are transgender (male to female) per km^2 PerTotalPop_Hetero 0.50 6.29E-149 8.93E-147 The sum of PLHIV who participate in heterosexual activities per 100,000 Sum_MF 0.49 5.18E-145 7.36E-143 The sum of all PLHIV who are transgender (male to female) Sum_Hisp 0.49 9.24E-144 1.31E-141 The sum of PLHIV who identify as hispanic ethnicity Top5_DR 0.48 1.02E-136 1.44E-134 The mean of the top five highest diversification rates PerTotalPop_FirstNation 0.48 5.98E-134 8.49E-132 The sum of PLHIV who identify as First Nations per 100,000 PerArea_MSM 0.47 4.28E-128 6.08E-126 The sum of PLHIV who identify as MSM per km^2 Sum_Female 0.47 1.04E-127 1.47E-125 The sum of all PLHIV who are female Top5_changeDR 0.46 2.09E-127 2.97E-125 The mean of the top five highest annual changes in diversification rates Sum_Black 0.46 1.25E-126 1.78E-124 The sum of PLHIV who identify as black ethnicity PerTotalPop_Black 0.46 2.32E-124 3.30E-122 The sum of PLHIV who identify as black ethnicity per 100,000 PerTotalPop_changeDR 0.44 2.18E-112 3.09E-110 The sum of all annual changes in diversification rate per 100,000 PerTotalPop_Hisp 0.44 1.43E-111 2.02E-109 The sum of PLHIV who identify as hispanic ethnicity per 100,000 PerArea_Hisp 0.43 1.54E-109 2.18E-107 The sum of PLHIV who identify as hispanic ethnicity per km^2 Top3_changeDR 0.43 2.60E-105 3.70E-103 The mean of the top three highest annual changes in diversification rates Sum_VL 0.42 2.70E-103 3.83E-101 The sum of all viral loads Top3_DR 0.42 3.69E-102 5.24E-100 The mean of the top three highest diversification rates PerTotalPop_Asian 0.42 6.01E-100 8.53E-98 The sum of PLHIV who identify as asian ethnicity per 100,000 PerTotalPop_VL 0.41 4.94E-99 7.02E-97 The sum of all viral loads per 100,000 Sum_Asian 0.41 1.39E-98 1.97E-96 The sum of PLHIV who identify as asian ethnicity PerArea_Asian 0.41 7.86E-96 1.12E-93 The sum of PLHIV who identify as asian ethnicity per km^2 PerTotalPop_Female 0.40 2.87E-92 4.07E-90 The sum of all PLHIV who are female per 100,000 PerArea_changeDR 0.40 1.30E-90 1.84E-88 The sum of all annual changes in diversification rate per km^2 PerArea_VL 0.37 6.50E-78 9.22E-76 The sum of all viral loads per km^2 Top5_VL 0.34 2.74E-64 3.89E-62 The mean of the top five highest viral loads Top3_VL 0.33 3.60E-63 5.12E-61 The mean of the top three highest viral loads Top5_log10VL 0.31 3.88E-54 5.50E-52 The mean of the top five highest log10 (viral loads) PerArea_FM 0.31 2.80E-52 3.98E-50 The sum of all PLHIV who are transgender (female to male) per km^2 Top5_lnDR 0.30 1.05E-51 1.48E-49 The mean of the top five highest ln (diversification rates) 114 Top5_log10DR 0.30 1.05E-51 1.48E-49 The mean of the top five highest log10 (diversification rates) Max_changeDR 0.30 1.96E-50 2.78E-48 The maximum annual change in diversification rate Top3_lnDR 0.30 4.71E-49 6.68E-47 The mean of the top three highest ln (diversification rates) Top3_log10DR 0.30 4.71E-49 6.68E-47 The mean of the top three highest log10 (diversification rates) PerPLHIVPop_NextFirstVL 0.29 2.45E-48 3.47E-46 The sum of new HIV cases in the subsequent year per 100 PLHIV Top3_log10VL 0.29 2.47E-46 3.50E-44 The mean of the top three highest log10 (viral loads) Max_VL 0.28 6.37E-44 9.05E-42 The maximum viral load Max_DR 0.27 3.50E-41 4.97E-39 The maximum diversification rate Max_log10DR 0.25 1.13E-34 1.61E-32 The maximum log10(diversification rate) Max_lnDR 0.25 1.13E-34 1.61E-32 The maximum ln(diversification rate) PerTotalPop_FM 0.25 4.03E-34 5.73E-32 The sum of all PLHIV who are transgender (female to male) per 100,000 Max_log10VL 0.25 9.00E-34 1.28E-31 The maximum log10(viral load) Sum_FM 0.22 1.73E-26 2.46E-24 The sum of all PLHIV who are transgender (female to male) PerPLHIVPop_MSM 0.10 3.98E-07 5.66E-05 The percent of PLHIV who identify as MSM PerPLHIVPop_Male 0.10 1.14E-06 1.62E-04 The percent of PLHIV who are male Median_log10VL 0.09 1.16E-05 1.64E-03 The median log10(viral load) PerPLHIVPop_FirstNation 0.08 4.00E-05 5.67E-03 The percent of PLHIV who identify as First Nations PerPLHIVPop_MF 0.08 1.21E-04 1.71E-02 The percent of PLHIV who are transgender (male to female) Median_log10DR 0.07 3.27E-04 4.65E-02 The median log10(diversification rate) Median_lnDR 0.07 3.27E-04 4.65E-02 The median ln(diversification rate) PerPLHIVPop_Hisp 0.07 8.63E-04 1.23E-01 The percent of PLHIV who identify as hispanice ethnicity Mean_log10VL 0.06 2.57E-03 3.65E-01 The mean log 10(viral load) Mean_changeDR 0.05 1.79E-02 1.00E+00 The mean change in diversification rate PerPLHIVPop_changeDR 0.05 2.30E-02 1.00E+00 The sum of all annual changes in diversification rate per 100 PLHIV PerPLHIVPop_log10VL 0.05 2.51E-02 1.00E+00 The sum of log10(viral load) per 100 PLHIV PerPLHIVPop_HCV 0.05 2.53E-02 1.00E+00 The percent of PLHIV who were ever coinfected with HCV PerPLHIVPop_BinaryChangeDR 0.04 6.22E-02 1.00E+00 The sum of PLHIV with an annual change in diversification rate >=1 per 100 PLHIV PerPLHIVPop_IDU 0.04 6.99E-02 1.00E+00 The percent of PLHIV who identify as IDU PerPLHIVPop_RiskKnown 0.03 1.34E-01 1.00E+00 The percent of PLHIV with some risk info known PerPLHIVPop_FirstVL 0.03 1.86E-01 1.00E+00 The sum of new HIV cases per 100 PLHIV PerPLHIVPop_White 0.02 2.88E-01 1.00E+00 The percent of PLHIV who identify as white ethnicity Mean_DR 0.01 4.80E-01 1.00E+00 The mean diversification rate PerPLHIVPop_DR 0.01 4.80E-01 1.00E+00 The sum of all diversification rates per 100 PLHIV PerPLHIVPop_Other 0.00 9.13E-01 1.00E+00 The percent of PLHIV who were exposed to other risks PerPLHIVPop_InfoKnown 0.00 9.73E-01 1.00E+00 The percent of PLHIV with all risk info known 115 PerPLHIVPop_Black 0.00 9.33E-01 1.00E+00 The percent of PLHIV who identify as black ethnicity PerArea_log10DR 0.00 8.37E-01 1.00E+00 The sum of all log10(diversification rates) per km^2 PerArea_lnDR 0.00 8.37E-01 1.00E+00 The sum of all ln(diversification rates) per km^2 Mean_Age -0.01 5.94E-01 1.00E+00 The mean age Mean_VL -0.01 5.79E-01 1.00E+00 The mean viral load Median_changeDR -0.01 5.13E-01 1.00E+00 The median annual change in diversification rate PerPLHIVPop_VL -0.01 4.95E-01 1.00E+00 The sum of all viral loads per 100 PLHIV PerPLHIVPop_Blood -0.02 4.62E-01 1.00E+00 The percent of PLHIV who were exposed to blood risk PerPLHIVPop_Asian -0.02 3.01E-01 1.00E+00 The percent of PLHIV who identify as asian ethnicity Median_DR -0.03 1.45E-01 1.00E+00 The median diversification rate PerPLHIVPop_FM -0.04 5.75E-02 1.00E+00 The percent of PLHIV who are transgender (female to male) Median_VL -0.05 1.80E-02 1.00E+00 The median viral load PerPLHIVPop_Hetero -0.07 9.25E-04 1.31E-01 The percent of PLHIV who participate in heterosexual activities Sum_lnDR -0.07 8.16E-04 1.16E-01 The sum of all ln(diversification rates) Sum_log10DR -0.07 8.16E-04 1.16E-01 The sum of all log10(diversification rates) PerPLHIVPop_Female -0.09 3.83E-06 5.44E-04 The percent of PLHIV who are female PerTotalPop_lnDR -0.10 1.21E-06 1.71E-04 The sum of all ln(diversification rates) per 100,000 PerTotalPop_log10DR -0.10 1.21E-06 1.71E-04 The sum of all log10(diversification rates) per 100,000 Mean_lnDR -0.19 2.49E-21 3.54E-19 The mean ln(diversification rate) Mean_log10DR -0.19 2.49E-21 3.54E-19 The mean log10(diversification rate) PerPLHIVPop_lnDR -0.19 2.49E-21 3.54E-19 The sum of all ln(diversification rates) per 100 PLHIV PerPLHIVPop_log10DR -0.19 2.49E-21 3.54E-19 The sum of all log10(diversification rates) per 100 PLHIV Table A-3 The goodness-of-fit of the full zero-inflated negative binomial (ZINB) model in comparison to competing models. The AIC, BIC, absolute log-likelihood, degrees of freedom of the ZINB model were compared to those of its zero-inflated Poisson (ZIP) and Hurdle equivalents, as well as a reduced ZINB with no diversification rate measures. Model AIC BIC |Log- likelihood| df Likelihood ratio test (relative to full ZINB), p-value Full ZINB 3,109 3,152 1,547 11 - Full ZIP 3,137 3,175 1,562 10 <0.001 Full Hurdle 3,128 3,170 1,556 11 <0.001 Reduced ZINB 3,238 3,265 1,614 8 <0.001 116 Table A-4 A comparison of the mean predictive power measures for the full zero-inflated negative binomial (ZINB) model to its reduced model without any diversification rate. The comparison criteria were calculated separately for each year in the training data (2009, 2010, 2012, 2013, 2015, and 2016) and testing data (2008, 2011, 2014) and the mean across those years is reported here. Model Pearson corr. (obs. cases >0) Sensitivity Specificity PPV NPV Accuracy Total pred. new cases Total obs. new cases TP FP TN FN Testing Dataset Full ZINB 0.77 0.66 0.68 0.61 0.73 0.67 176.3 218.0 70.3 45.7 97.7 36.3 Reduced ZINB 0.66 0.75 0.61 0.59 0.76 0.67 213.3 218.0 79.7 56.3 87.0 27.0 Training Dataset Full ZINB 0.67 0.72 0.63 0.57 0.77 0.66 191.5 204.8 72.5 55.5 93.3 28.7 Reduced ZINB 0.58 0.77 0.56 0.54 0.78 0.64 201.2 204.8 77.7 65.8 83.0 23.5 *Criteria include the Pearson’s correlation coefficient (for observed cases > 0 to remove the effect of zero-inflation); sensitivity (true positives/ (true positives + false negatives)); specificity (true negatives/(true negatives + false positives)); positive predictive value (true positives/ (true positives + false positives)); negative predictive value (true negatives/(true negatives + false negatives)); total number of predicted new cases; total number of observed new cases; true positives; false positives; true negatives; and false negatives. Table A-5 The annual number of new HIV cases with sequences compared to those with census tract (CT) information available and estimated provincial incidence. Data was last updated on February 4 2019. Year New HIV cases with sequences New HIV cases with sequences and CT info Percent new HIV cases with CT info BC Population New HIV cases with sequences in BC per 100,000 2008 379 239 63.06 4,349,412 8.71 2009 374 253 67.65 4,410,679 8.48 2010 364 226 62.09 4,465,924 8.15 2011 363 223 61.43 4,499,139 8.07 2012 313 202 64.54 4,546,290 6.88 2013 310 200 64.52 4,590,081 6.75 2014 293 210 71.67 4,646,462 6.31 2015 308 214 69.48 4,694,699 6.56 2016 294 213 72.45 4,757,658 6.18 2017 211 163 77.25 4,817,160 4.38 2018 201 64 31.84 4,849,442 4.14 117 A.3 Supplementary data and shapefiles Statistics Canada census boundary files Statistics Canada. 2016. lct_000a16a_e. 2016 Census. Statistics Canada Catalogue no. 92-168-XWE. Released November 16, 2016 https://www12.statcan.gc.ca/census-recensement/alternative_alternatif.cfm?l=eng&dispext =zip&teng=lct_000a16a_e.zip&k=%20%20%20%20%207190&loc=http://www12.statcan.gc.ca/census-recensement/2011/geo/bound-limit/files-fichiers/2016/lct_000a16a_e.zip (accessed August 9, 2017). Statistics Canada census tract data Statistics Canada. Release date: February 8, 2017. Census Profile, 2016 Census. Last updated May 30, 2018. Catalogue no. 98-316-X2016001 https://www12.statcan.gc.ca/census-recensement/2016/dp-pd/prof/index.cfm?Lang=E (accessed October 2, 2017). Statistics Canada. Release date: October 24, 2012. 2011 Census Profile, 2011 Census. Last updated: Date modified: 2018-05-23. Catalogue no. 98-316-XWE. https://www12.statcan.gc.ca/census-recensement/2011/dp-pd/prof/index.cfm?Lang=E (accessed October 2, 2017). Statistics Canada. Release date: March 13, 2007. 2006 Community Profiles, 2006 Census. Last updated: January 5, 2015. Catalogue no. 92-591-XWE. https://www12.statcan.gc.ca/census-recensement/2006/dp-pd/prof/92-591/index.cfm?Lang=E (accessed October 2, 2017). 118 Statistics Canada. Release date: 2002. Population and Dwelling Counts, for Canada, Provinces and Territories, 2001 and 1996 Censuses - 100% Data, 2001 Census. Last updated: 2012-07-04. Catalogue no. n/a. https://www12.statcan.gc.ca/english/census01/products/standard/popdwell/Table-PR.cfm (accessed October 2, 2017). Statistics Canada postal code conversion file 2018-10-10, "Postal Code Conversion File, August 2018 Postal Codes, 2018", http://hdl.handle.net.ezproxy.library.ubc.ca/11272/10659 Statistics Canada [Distributor] V2 [Version]. Health authority boundary files BC Stats. Release date: December 11, 2014. Health Authority Boundaries. Last updated: November 3, 2017. https://catalogue.data.gov.bc.ca/dataset/health-authority-boundaries (accessed February 15, 2019).
Thesis/Dissertation
eng
Population and Public Health
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Applying phylogenetic lineage-level diversification rate to infer epidemiological risk factors and spatiotemporal distribution of HIV transmission in British Columbia
Text