Www.medrxiv.org



Table A1: Parameter definitions and sources.ParameterDefinitionData sourceValuebetaBasic transmission rate SPIM psi.spim <- 0.2 # epidemic growth rate (from PB 26/03) dI.spim <- 4.6 # mean infectious period (from PB 26/03) dL.spim <- 2 # mean latent period (from PB 26/03) R0.spim <- 2.5 #(from PB 26/03) beta <- R0.spim/(dL.spim + dI.spim) # transmission rate per contacttrTesting rate of symptomatic in hospitalAssumed all symptomatics will be tested within 2 daysn (time taken for test result to come back)iRate of transition from Exposed (asymptomatic) -> I (symptomatic)SPIM2 daysrsevereRecovery rate of those infected with COVID19 in hospital (regardless of whether they picked it up in the community or nosocomially) , , - 13.5 dayshsRecruitment rate of susceptible patients1 - hi - he - hrheRecruitment rate of asymtpomatic patientsR0*hihrRecruitment rate of recovered patients0.00Assume no recruitment of recovered individualsatProportion of infected patients picked up at admission COVID – Approaches to Reducing Nosocomial Infection Questionnaire Cleaned v2.xlsx, Sheet 2 (positive_tests_0904.xlsx) has cumulative number of inpatients tested from 48hrs after admission (these are in patient tests) on 09/04. Total infected admissions - those tested after 48hrs = those that were infected but missed0.8atrRandom admissions testing rateCOVID – Approaches to Reducing Nosocomial Infection Questionnaire Cleaned v2.xlsx, Sheet 2 (positive_tests_0904.xlsx) - Ratio of positive tests at admission to total tests done, scaled for daily imported cases 0.10?dsDischarge rate susceptible0.10deDischarge rate exposed0.05diDischarge rate infected0.05drDischarge rate recovered1/time taken for negative test result0.50msDeath rate susceptibleNHS average in hospital death rate0.0045meDeath rate exposed average rate/average time infected (from PHE tracker)0.064/rsevere miDeath rate infected average rate/average time infected (from PHE tracker)0.064/rsevere mrDeath rate recoveredmseRelative infeciousness of exposed vs infected patient1.00nRate of negative test result coming back0.50cRelative infeciousness of cohort isolation1.00Assumption: isolation is 100% effective and that staff pick up infections from uncohorted/isolated patients. When we have data about infections per ward we can change this and recalibrate but for basic UK calibration it doesn’t matter much as long as the total number of infections are correct. aabsence proportion of HCWs ADDIN ZOTERO_ITEM CSL_CITATION {"citationID":"AyKqfMeM","properties":{"formattedCitation":"(Kluytmans et al., 2020)","plainCitation":"(Kluytmans et al., 2020)","noteIndex":0},"citationItems":[{"id":449,"uris":[""],"uri":[""],"itemData":{"id":449,"type":"article-journal","title":"SARS-CoV-2 infection in 86 healthcare workers in two Dutch hospitals in March 2020","container-title":"medRxiv","page":"2020.03.23.20041913","source":"","abstract":"<p><b>Background</b> On February 27, 2020, the first patient with COVID-19 was reported in the Netherlands. During the following weeks, nine healthcare workers (HCWs) were diagnosed with COVID-19 in two Dutch teaching hospitals, eight of whom had no history of travel to China or Northern-Italy. A low-threshold screening regimen was implemented to determine the prevalence and clinical presentation of COVID-19 among HCWs in these two hospitals. <b>Methods</b> HCWs who suffered from fever or respiratory symptoms were voluntarily tested for SARS-CoV-2 by real-time reverse-transcriptase PCR on oropharyngeal samples. Structured interviews were conducted to document symptoms for all HCWs with confirmed COVID-19. <b>Findings</b> Thirteen-hundred fifty-three (14%) of 9,705 HCWs employed were tested, 86 (6%) of whom were infected with SARS-CoV-2. Most HCWs suffered from relatively mild disease and only 46 (53%) reported fever. Seventy-nine (92%) HCWs met a case definition of fever and/or coughing and/or shortness of breath. None of the HCWs identified through the screening reported a travel history to China or Northern Italy, and 3 (3%) reported to have been exposed to an inpatient known with COVID-19 prior to the onset of symptoms. <b>Interpretation</b> Within two weeks after the first Dutch case was detected, a substantial proportion of HCWs with fever or respiratory symptoms were infected with SARS-CoV-2, probably caused by acquisition of the virus in the community during the early phase of local spread. The high prevalence of mild clinical presentations, frequently not including fever, asks for less stringent use of the currently recommended case-definition for suspected COVID-19.</p>","DOI":"10.1101/2020.03.23.20041913","language":"en","author":[{"family":"Kluytmans","given":"Marjolein"},{"family":"Buiting","given":"Anton"},{"family":"Pas","given":"Suzan"},{"family":"Bentvelsen","given":"Robbert"},{"family":"Bijllaardt","given":"Wouter","dropping-particle":"van den"},{"family":"Oudheusden","given":"Anne","dropping-particle":"van"},{"family":"Rijen","given":"Miranda","dropping-particle":"van"},{"family":"Verweij","given":"Jaco"},{"family":"Koopmans","given":"Marion"},{"family":"Kluytmans","given":"Jan"}],"issued":{"date-parts":[["2020",3,31]]}}}],"schema":""} Kluytmans et al., (2020)37%tHCWtesting rate of HCWs0No testing of asymptomatic HCWs in regular scenariobP2PRelative infectiousness P2PThe scaling factor of the transmission rate, β for spread between patient to patient, patient to HCW, HCW to HCW, and HCW to patient, and hi, the proportion of daily cases admitted to the hospital were calibrated to reflect the average proportion of HCWs that tested positive from personal correspondence from two English Trusts0.25bP2HRelative infectiousness P2H0.045bH2PRelative infectiousness H2P0.5bH2HRelative infectiousness H2H0.4hiRecruitment rate of symtpomatic patientsNHS situation report daily COVID19 admissions0.05z_capPositive cohort capacityAssumed unlimited?mzDeath rate positive cohortmiprevalenceInitial population prevalencePHE dashboard datazRelative infectiousness within isolation (to HCWs)yRelative infectiousness within confirmed negative cohort doubling_timeDoubling time of epidemic SPI-M report2.6rhealthRecovery rate of HCWs infected with COVID19SPI-M report1/irabsentRecover rate of absent HCWs7 days Supplementary File 1Estimating the average number of SARS-CoV2 infections in different types of people in hospital settings using the approach of Heesterbeek and Roberts (2010)The approach of Heesterbeek and Roberts (2010) was used to calculate the Next-Generation Matrix (NGM) for transmission between hospitalised people E, those suspected to have infection (Ex), those who have tested negative (Ey) and health care workers (Ehcw). Using the notation by Heesterbeek, the row labels reflect the recipient of the infection and the column label refers to the person doing the infecting:NGM, K, assuming no contact tracing, is;EExEyEhcwE0.2420.0000.0006.91 × 10?5Ex0.0680.1700.0686.91 × 10?5Ey0.0000.0000.1213.45 × 10?5Ehcw1.38 × 10?43.27 × 10?41.35 × 10?40.452The within-hospital R0 is 0.452 and is represented by the dominant eigenvalue of the NGM. The number of secondary infections from each source are as follows:From a HCW to any type of patient = 1.73 × 10?4 = 6.91 × 10?5 + 6.91 × 10?5 + 3.45 × 10?5From any patient to any type of patient = 0.31 (= 0.242 + 0.068), assuming that suspect patients make up a small proportion of the population. This value goes down to 0.17 if the proportion of hospitalised patients that are suspect increases.From any patient to any HCW = 1.38 × 10?4, assuming that suspect and test-negative patients make up a small proportion of the population. This increases to 0.008 if, pessimistically, all patients were to be suspect patients 3.27 × 10?4.From a typical HCW to a HCW is 0.452.Among suspect patients 0.17.From test-negative to suspect patients, 0.068.Among test-negative patients, 0.12.The remaining pages show the matrices produced when generating the NGM, and the full R code follows: Matrix T: EIExIxEyIyZEhcwIhcw9.47 × 10?29.47 × 10?20.000.000.000.000.002.32 × 10?52.32 × 10?50.000.000.000.000.000.000.000.000.000.000.009.47 × 10?29.47 × 10?20.000.000.002.32 × 10?52.32 × 10?50.000.000.000.000.000.000.000.000.000.000.000.000.004.73 × 10?24.73 × 10?20.001.16 × 10?51.16 × 10?50.000.000.000.000.000.000.000.000.000.000.000.000.000.000.000.000.000.001.70 × 10?61.70 × 10?61.14 × 10?41.14 × 10?48.51 × 10?78.51 × 10?71.70 × 10?67.58 × 10?27.58 × 10?20.000.000.000.000.000.000.000.000.00Matrix Sigma: EIExIxEyIyZEhcwIhcw?5.75 × 10?10.000.000.000.000.000.000.000.003.03 × 10?1?6.47 × 10?10.000.000.000.000.000.000.000.000.00?8.58 × 10?10.000.000.000.000.000.000.005.00 × 10?13.03 × 10?1?5.64 × 10?10.005.00 × 10?10.000.000.000.000.000.000.00?5.75 × 10?10.000.000.000.000.000.000.000.003.03 × 10?1?6.47 × 10?10.000.000.000.000.005.00 × 10?15.00 × 10?10.000.00?1.05 × 10?10.000.000.000.000.000.000.000.000.00?5.21 × 10?10.000.000.000.000.000.000.000.001.91 × 10?1?3.45 × 10?1Matrix Sigma Inverse:EIExIxEyIyZEhcwIhcw1.74?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.008.15 × 10?11.55?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.001.17?0.00?0.00?0.00?0.00?0.00?0.007.23 × 10?11.376.27 × 10?11.777.23 × 10?11.37?0.00?0.00?0.00?0.00?0.00?0.00?0.001.74?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.008.15 × 10?11.55?0.00?0.00?0.003.436.508.508.413.436.509.49?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.001.92?0.00?0.00?0.00?0.00?0.00?0.00?0.00?0.001.062.90Matrix KL (NGM large domain):EIExIxEyIyZEhcwIhcw2.42 × 10?11.46 × 10?10.000.000.000.000.006.91 × 10?56.71 × 10?50.000.000.000.000.000.000.000.000.006.84 × 10?21.30 × 10?11.70 × 10?11.68 × 10?16.84 × 10?21.30 × 10?10.006.91 × 10?56.71 × 10?50.000.000.000.000.000.000.000.000.000.000.000.000.001.21 × 10?17.32 × 10?20.003.45 × 10?53.35 × 10?50.000.000.000.000.000.000.000.000.000.000.000.000.000.000.000.000.000.009.23 × 10?51.69 × 10?42.18 × 10?42.16 × 10?49.01 × 10?51.68 × 10?41.61 × 10?52.26 × 10?12.19 × 10?10.000.000.000.000.000.000.000.000.00Auxiliary Matrix E:EExEyEhcw1.000.000.000.000.000.000.000.000.001.000.000.000.000.000.000.000.000.001.000.000.000.000.000.000.000.000.000.000.000.000.001.000.000.000.000.00R script:# Estimating the number of secondary infections in different types of people using the # approach of Heesterbeek and Roberts. The parameters until lines 170 are taken directly from the code load_parameters_and_ics_3cohort_transm_uk.R# The matrix calculations are from line 170 onwards. # Corrected 26.04.20 to include a missing term in matrix Sigma# Parameters taken from file load_parameters_and_ics_3cohort_transm_uk.R#New parameters:y <- 0.5 #Sy.t0 -see below#Ey.t0 -see below#Iy.t0 -see below#Ry.t0 -see below#tot_y.t0 -see belowN <- 800000/2hosp_cap <- 668z_cap <- hosp_cap # unlimited capacity#hosp_cap = hosp_cap - z_cap## spi-m-o report parameterspsi.spim <- 0.2 # epidemic growth rate (from PB 26/03)dI.spim <- 4.6 # mean infectious period (from PB 26/03)dL.spim <- 2 # mean latent period (from PB 26/03)R0.spim <- 2.5 #(from PB 26/03)beta <- R0.spim/(dL.spim + dI.spim) # transmission rate per contactL <- dL.spimi <- exp(-1/2)*1/2 # rate of infection (1/latent period)rsevere <- exp(-1/11)*1/11 # rate of recoveryrhealth <- 1/dI.spimrabsent <- 1/7## if this is changed it needs to move into the parameter function# iso_efficacy = sample(1:4, 1)iso_efficacy = -1if(iso_efficacy == -1){ z = 1 ## from hci paper c = 1 } else if(iso_efficacy ==1){ z = 0 # both 100 % effective c = 0} else if(iso_efficacy ==2){ z = 0 c = 0.5} else if(iso_efficacy ==3){ z = 1 # not effective at all c = 1 # not effective at all} else if(iso_efficacy ==4){ z = 0.5 c = 0.75}prevalence <- 1.2*10^-6 # low and high prevalencedoubling_time <- 2.6hi <- 0.05 # this is multiplied by prevalence in the model_function script so is proportion of community cases assigned to this hosptial# capacity/population (trust data + fingertips)he <- R0.spim*hi # assume R0*hi people are exposed.# hr <- sample(c(0.2, 0.5, 0.8),1) - he # assume different proportions of the population have already had ithr <- 0hs <- 1 - hr - heif(hs < 0){ print("hs less than zero!!!!")}I.t0<-0 # infectedE.t0<-0# exposedR.t0<-0 # recoveredS.t0 <- hosp_cap - E.t0 - I.t0 - R.t0if(hi*(hosp_cap + z_cap)< z_cap){ isolated.t0 = hi*(hosp_cap + z_cap) cohort.t0 = 0} else { isolated.t0 = z_cap cohort.t0 = hi*(hosp_cap + z_cap) - z_cap}## at least one of these has to be bigger than 0Sx.t0<- 10 # Susceptible in suspected cohortEx.t0<-0 # Exposed in suspected cohortIx.t0<-0# Infected in suspected cohortRx.t0<-0 # Recovered in suspected cohortZ.t0<- 0 # tested +ive cohortSy.t0<-10 # Susceptible in tested -ive cohortEy.t0<-0 # Exposed in tested -ive cohortIy.t0<-0# Infected in tested -ive cohortRy.t0<-0 # Recovered in tested -ive cohortn_staff = 8180# initial values - HCWsEhcw.t0<- R0.spim*prevalence*n_staff # exposedIhcw.t0<-prevalence*n_staff # infectedShcw.t0<-n_staff - Ehcw.t0 # susceptible (4 patients per staff member)?)Rhcw.t0<-0 # recoveredAhcw.t0<-0 # absenttot_hcw.t0<-sum(c(Shcw.t0,Ehcw.t0,Ihcw.t0,Rhcw.t0))S.t0 <- hosp_cap - E.t0 - I.t0 - R.t0 -Sx.t0 - Ex.t0 - Ix.t0 - Rx.t0 - Sy.t0 - Ey.t0 - Iy.t0 - Ry.t0 - Z.t0 # susceptibleP.t0 <- N - S.t0 - E.t0 - I.t0 - R.t0 - Sx.t0 - Ex.t0 - Ix.t0 - Rx.t0 - Sy.t0 - Ey.t0 - Iy.t0 - Ry.t0 - Z.t0 # communitytot.t0<-sum(c(S.t0,E.t0,I.t0,R.t0,Sx.t0,Ex.t0,Ix.t0,Rx.t0,Sy.t0,Ey.t0,Iy.t0,Ry.t0,P.t0,Z.t0))tot_hosp.t0<-sum(c(S.t0,E.t0,I.t0,R.t0,Sx.t0,Ex.t0,Ix.t0,Rx.t0,Sy.t0,Ey.t0,Iy.t0,Ry.t0,Z.t0))tot_x.t0=sum(c(Sx.t0,Ex.t0,Ix.t0,Rx.t0))tot_y.t0=sum(c(Sy.t0,Ey.t0,Iy.t0,Ry.t0))# General parameters# testingat <- 0.8 # proportion of infected (symptomatic) individuals tested are tested at admission (and therefore get +ve results and isolated)atr <- 0.1 # proportion of admissions tested at random (influenza etc).. taken from admissions data during winter for RTI and flu, might not be accurate nowatrE <- 0.1 # proportion of exposed hospitalisations testing on admission (same as 'atr' for now)ct <- 0ctHCW <- 0# unknown/estimatedds <- 1/10 # rate of discharge de <- 0.5*ds # rate of discharge (from exposed comp.)di <- 0.5*ds # rate of discharge (from infected comp.)dr <- 0.5 # rate of discharge (from recovered comp.) go home when they are recoveredms <- 0.0045 # mortality rate (from susceptible comp.) (NHS average in hospital death rate)mi <- 0.064/(dI.spim) # mortality rate (from infected comp.) average rate/average time infected (from PHE tracker)me <- ms # mortality rate (from exposed comp.)mz <- mimr <- ms # mortality rate (from recovered comp.)e <- 1 # relative infectiousness of exposed vs. infected (0<=e<=1)n <- 0.5 # rate of negative test result (1/time to result)# how long does it take for in-hospital symptomatics to get diagnosed?tr <- n # rate of testing (of people that are left). test everyone on the ward that are infected within 2 days# -hcw parameters# red_bP2P <- sample(c(0,0.5,0.75,1), 1)a <- 0.37red_bP2P <- 0red_bP2H <- 0red_bH2P <- 0red_bH2H <- 0tHCW <- n*exp(-(1- 0.63))*(1-0.63) # testing rate for infected HCW (subsequenly isolated.. 2 days)# rate in dutch hospital preprint 86/9705 hcws over 10 days with 5*15 contacts per day assuming staff only work 5/7 days. / beta/number of contacts in general pop since beta is infection rate*contacts in general popbP2P <- 0.25*(1 - red_bP2P) # relative infection rate patient to patient (assume 4 contacts a day)# bP2H <- inf_rate*(1 - red_bP2H) # control of relative infectiousness for patient to HCW# bH2P <-0.05*inf_rate*(1 - red_bH2P) # control of relative infectiousness for HCW to patientbH2H <- 0.2 # control of relative infectiousness for HCW to HCWbP2H <- 0.003 # control of relative infectiousness for patient to HCWbH2P <- 0.5 # control of relative infectiousness for HCW to patient# period_to_test <- sample(c(1,7,14, 100), 1)period_to_test = -1ics <-c(S.t0, E.t0, I.t0*runif(1, 0.5,1.5), R.t0, Sx.t0, Ex.t0, Ix.t0, Rx.t0,Sy.t0, Ey.t0, Iy.t0, Ry.t0, Z.t0, P.t0, tot.t0, tot_hosp.t0, Shcw.t0, Ehcw.t0, Ihcw.t0, Rhcw.t0, Ahcw.t0, tot_hcw.t0, tot_x.t0, tot_y.t0)params <- c(beta, tr, ct, i, rsevere, hs, he*runif(1, 0.5,1.5), hi, hr, at, atr, atrE, ds, de, di, dr*runif(1, 0.5,1.5), ms, me, mi, mr, e, n, c, a*runif(1, 0.5,1.5), tHCW, ctHCW, bP2P, bP2H, bH2P, bH2H, z_cap, mz, rgamma(1, 1.43, 55768), z, red_bH2P, red_bH2H, red_bP2H, red_bP2P, doubling_time, rhealth,y*runif(1,0.5,1.5), runif(1,7,14), period_to_test) names(ics) <- c("S.t0","E.t0","I.t0","R.t0","Sx.t0","Ex.t0","Ix.t0","Rx.t0","Sy.t0","Ey.t0","Iy.t0","Ry.t0","Z.t0","P.t0","tot.t0","tot_hosp.t0","Shcw.t0","Ehcw.t0","Ihcw.t0","Rhcw.t0","Ahcw.t0","tot_hcw.t0","tot_x.t0","tot_y.t0")names(params) <- c("beta","tr","ct","i","r","hs","he","hi","hr","at","atr","atrE","ds","de","di","dr","ms","me","mi","mr","e","n","c","a","tHCW","ctHCW","bP2P","bP2H","bH2P","bH2H","z_cap","mz","prevalence","z","red_bH2P","red_bH2H","red_bP2H","red_bP2P","doubling_time","rhealth","y", "rabsent", "period_to_test") # Set up matrix T, following Diekmann, Heesterbeek and Roberte (2010) J Roy SOc Interfaces# For now, taking totals in the suspected, test negative and positive people to equal the # values for Sx.t0 and Sy.t0, as specified above. Can set these to zero to turn off effect # of cohorting (though the corresponding transmission rates will need to be dropped if they're being# divided by these values)# or just ignore the values relating to these people in the NGMNx <- Sx.t0 # total in suspectedNy <- Sy.t0 # total in test-negativeTrEfromE <- beta*e*bP2PTrEfromI <- beta* bP2P(TrEfromEhcw <- beta*e*bH2P/n_staff)(TrEfromIhcw <- beta* bH2P/n_staff)(TrExfromEx <- beta*c*e*bP2P)(TrExfromIx <- beta*c*bP2P)(TrExfromEhcw <- beta*c*e*bH2P/n_staff)(TrExfromIhcw <- beta*c* bH2P/n_staff)(TrEyfromEy <- beta*y* e*bP2P)(TrEyfromIy <- beta*y* bP2P)(TrEyfromEhcw <- beta*y* e*bH2P/n_staff)(TrEyfromIhcw <- beta*y* bH2P/n_staff)(TrEhcwfromE <- beta*e*bP2H/hosp_cap)(TrEhcwfromI <- beta*bP2H/hosp_cap)(TrEhcwfromEx <- beta* c* e*bP2H/Nx)(TrEhcwfromIx <- beta* c*bP2H/Nx)(TrEhcwfromEy <- beta*y*e*bP2H/hosp_cap )(TrEhcwfromIy <- beta*y*bP2H/hosp_cap)(TrEhcwfromZ <- beta*z*bP2H/hosp_cap)(TrEhcwfromEhcw <- beta*e*bH2H)(TrEhcwfromIhcw <- beta*bH2H)cyz#T <- matrix( c(TrEfromE, TrEfromI,0,0,0,0,0,TrEfromEhcw,TrEfromIhcw, 0,0,0,0,0,0,0,0,0, 0,0,TrExfromEx,TrExfromIx,0,0,0,TrExfromEhcw,TrExfromIhcw, 0,0,0,0,0,0,0,0,0, 0,0,0,0,TrEyfromEy,TrEyfromIy,0,TrEyfromEhcw,TrEyfromIhcw, 0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0, TrEhcwfromE,TrEhcwfromI,TrEhcwfromEx,TrEhcwfromIx, TrEhcwfromEy,TrEhcwfromIy,TrEhcwfromZ,TrEhcwfromEhcw,TrEhcwfromIhcw, 0,0,0,0,0,0,0,0,0), nrow=9,ncol=9, byrow=TRUE)dimnames(T)= list(c("E","I","Ex","Ix","Ey","Iy","Z","Ehcw","Ihcw"), c("E","I","Ex","Ix","Ey","Iy","Z","Ehcw","Ihcw"))T# Set up matrix Sigma for the transition rates (movements) between states, following # Diekmann, Heesterbeek and Roberte (2010) J Roy SOc Interfaces# In the following mXY holds the coefficient for Y that's used in the differential equation for X# So mEE holds the coefficient in from of E that's used in the differential equation for E# For now, setting the contact testing rate, ctr, to equal zero. In the model, it's calculated as # n*infected/susceptible, where n is the number of contacts tested. The same applies for ctrHCW.ctr <- 0ctrHCW <-0(mEE <- - me - de - ctr - i - rhealth)(mIE <- i)(mII <- - mi - rsevere - tr - di)(mExE <- ctr)(mExEx <- - me - i - de - n)(mExEy <- ctr)(mIxI <- tr)(mIxEx <- i )(mIxIx <- - mi - di - n )(mIxIy <- tr)(mEyEy <- - me - i - de - rhealth - ctr)(mIyEy <- i)(mIyIy <- - mi- di - tr - rsevere)(mZEx <- n)(mZIx <- n)(mZZ <- - mz - di - n*rsevere )(mEhcwEhcw <- - i - ctrHCW - rhealth)(mIhcwEhcw <- i*(1-a))(mIhcwIhcw <- - rhealth - tHCW)Sigma <- matrix( c(mEE,0,0,0,0,0,0,0,0, mIE,mII,0,0,0,0,0,0,0, mExE,0,mExEx,0,mExEy,0,0,0,0, 0,mIxI,mIxEx,mIxIx,0,mIxIy,0,0,0, 0,0,0,0,mEyEy,0,0,0,0, 0,0,0,0,mIyEy,mIyIy,0,0,0, 0,0,mZEx,mZIx,0,0,mZZ,0,0, 0,0,0,0,0,0,0,mEhcwEhcw,0, 0,0,0,0,0,0,0,mIhcwEhcw,mIhcwIhcw), nrow=9,ncol=9, byrow=TRUE)dimnames(Sigma)= list(c("E","I","Ex","Ix","Ey","Iy","Z","Ehcw","Ihcw"), c("E","I","Ex","Ix","Ey","Iy","Z","Ehcw","Ihcw"))Sigma(minusSigmaI <- -solve(Sigma))# Calculate the NGM with large domain K_L following # Diekmann, Heesterbeek and Roberte (2010) J Roy Soc Interfaces. This is calculated as# K_L = -T * SigmaI, or equivalently, T*minusSigmaI(K_L <- T %*% minusSigmaI)# Following Diekmann, Heesterbeek and Roberte (2010) J Roy Soc Interfaces, calculate the NGM with no non-zero # rows, K, using the equation K= E'*K_L*E, where E consists of the unit column vectors, e_i, for which the ith row of# T is not zero. In our case, the 1st, 3rd, 5th, 8th rows are non-zero, so E needs to be a 9x4 matrix with 1 in the # 1st, 3rd, 5th and 8th rows and zeros elsewhere. E' denotes the transpose of E.E <- matrix(c(1,0,0,0, 0,0,0,0, 0,1,0,0, 0,0,0,0, 0,0,1,0, 0,0,0,0, 0,0,0,0, 0,0,0,1, 0,0,0,0), nrow=9,ncol=4,byrow=TRUE)dimnames(E) <- list(c("E","I","Ex","Ix","Ey","Iy","Z","Ehcw","Ihcw"), c("E","Ex","Ey","Ehcw"))EEtransp <- t(E)(K <- Etransp %*% K_L %*% E) # NGM(R0 <- max(eigen(K)$values))library(dplyr)library(tibble)dfk = as.data.frame(K)rowSums(dfk)gt(as.data.frame(dfk) %>% rownames_to_column(" "))tab = gt(as.data.frame(dfk)%>% rownames_to_column(" ")) %>% fmt_scientific(columns = colnames(dfk), decimals = 2) %>% fmt_number(colnames(dfk)[1:3], rows = 1:3, decimals = 3) Model Functions######################################## COVID19 transmission model in hospital with 3 cohorts - suspcted +ive (X), tested -ive (Y) and tested +ive (Z)# with transmission in X and Y########################################rm(list=ls())require("deSolve")require("manipulate")# **********************************************************************************## FUNCTION: mod.dynmod.dyn <- function(t,var,par) { # Basic model with transmission in the hospital # And movement of patients from the community # Suspected +ive are held in X awaiting test results, -ives moved to Y, and +ives moved to Z (uncapped) # Transmission in cohorts X and Y # S<- var[["S.t0"]] # Susceptible to colonisation E<-var[["E.t0"]] # Exposed to COVID - asymptomatic I<-var[["I.t0"]] # Infected with COVID R<-var[["R.t0"]] # Recovered Sx<-var[["Sx.t0"]] # Susceptible in suspected cohort Ex<-var[["Ex.t0"]] # Exposed in suspected cohort Ix<-var[["Ix.t0"]] # Infected in suspected cohort Rx<-var[["Rx.t0"]] # Recovered in suspected cohort Sy<-var[["Sy.t0"]] # Susceptible in tested -ive cohort Ey<-var[["Ey.t0"]] # Exposed in tested -ive cohort Iy<-var[["Iy.t0"]] # Infected in tested -ive cohort Ry<-var[["Ry.t0"]] # Recovered in tested -ive cohort Z<-var[["Z.t0"]] # tested +ive cohort P<-var[["P.t0"]] # Community - pool entering/leaving hospital tot<-var[["tot.t0"]] # Total population tot_hosp<-var[["tot_hosp.t0"]] # Total hospitalised population Shcw<-var[["Shcw.t0"]] # Susceptible HCW Ehcw<-var[["Ehcw.t0"]] # Exposed HCW Ihcw<-var[["Ihcw.t0"]] # Infected HCW Rhcw<-var[["Rhcw.t0"]] # Recovered HCW Ahcw<-var[["Ahcw.t0"]] # Absent HCW tot_hcw<-var[["tot_hcw.t0"]] # total available HCW tot_x<-var[["tot_x.t0"]] # total patients in suspected cohort tot_y<-var[["tot_y.t0"]] # total patients in tested -ive cohort # model parameters for hospital spread beta<-par[["beta"]] # transmission rate tr<-par[["tr"]] # rate of testing ct<-par[["ct"]] # rate of contact testing i<-par[["i"]] # rate of infection (1/latent period) rsevere<-par[["r"]] # rate of recovery rhealth <-par[["rhealth"]] # recovery of HCW rabsent <-par[["rabsent"]] # return to work of HCW hs<-par[["hs"]] # proportion of hospitalisation into susceptible comp. (1=hs+he+hi+hr) he<-par[["he"]] # proportion of hospitalisation into exposed comp. (1=hs+he+hi+hr) hi<-par[["hi"]] # proportion of hospitalisation into infected comp. (1=hs+he+hi+hr) hr<-par[["hr"]] # proportion of hospitalisation into recovered comp. (1=hs+he+hi+hr) at<-par[["at"]] # proportion of infected hospitalisations that are tested on admission atr<-par[["atr"]] # proportion of non-infectious hospitalisations testing on admission atrE<-par[["atrE"]] # proportion of exposed hospitalisations testing on admission ds<-par[["ds"]] # rate of discharge (from susceptible comp.) de<-par[["de"]] # rate of discharge (from exposed comp.) di<-par[["di"]] # rate of discharge (from infected comp.) dr<-par[["dr"]] # rate of discharge (from recovered comp.) ms<-par[["ms"]] # mortality rate (from susceptible comp.) me<-par[["me"]] # mortality rate (from exposed comp.) mi<-par[["mi"]] # mortality rate (from infected comp.) mr<-par[["mr"]] # mortality rate (from recovered comp.) mz<-par[["mz"]] e<-par[["e"]] # relative infectiousness of exposed vs. infected (0<=e<=1) n<-par[["n"]] # rate of receiving test result (1/time to getting result) c<-par[["c"]] # relative infectiousness of X cohorting vs. infected (0<=c<=1) y<-par[["y"]] # relative infectiousness of Y cohorting vs. infected (0<=y<=1) z<-par[["z"]] # relative infectiousness of Z cohorting vs. infected (0<=z<=1) a<-par[["a"]] # absence proportion of all HCW infected with COVID-19 tHCW<-par[["tHCW"]] # testing rate for HCW ctHCW<-par[["ctHCW"]] # contact testing rate for HCW bP2P<-par[["bP2P"]] # control of relative infectiousness for patient to patient bP2H<-par[["bP2H"]] # control of relative infectiousness for patient to HCW bH2P<-par[["bH2P"]] # control of relative infectiousness for HCW to patient bH2H<-par[["bH2H"]] # control of relative infectiousness for HCW to HCW z_cap<-par[["z_cap"]] population_prevalence<-par[["prevalence"]] red_bH2P<-par[["red_bH2P"]] red_bH2H<-par[["red_bH2H"]] red_bP2H<-par[["red_bP2H"]] red_bP2P<-par[["red_bP2P"]] doubling_time<-par[["doubling_time"]] rhealth <-par[["rhealth"]] # recovery of HCW period_to_test <- par[["period_to_test"]] #double prevalence every 3.4 days prop_admitted <- hi*calculate_incidence(t - 2, population_prevalence) population_prevalence <- calculate_incidence(t, population_prevalence) tHCW = testing_rate_hcw(round(t), period_to_test, tHCW) a = absent_rate_hcw(round(t), period_to_test, a) N <- S+E+I+R+Sx+Ex+Ix+Rx+Sy+Ey+Iy+Ry+Z+P # total population size Nh <- S+E+I+R+Sx+Ex+Ix+Rx+Sy+Ey+Iy+Ry+Z# total hospital population size Nhcw <- Shcw+Ehcw+Ihcw+Rhcw # total available HCW population Nx <- Sx+Ex+Ix+Rx # total patients in suspected cohort Ny <- Sy+Ey+Iy+Ry # total patients in tested -ive cohort # Constant popn size in hospital discharge <- ds*S + de*E + di*I + dr*R + di*Z + di*Ix + ds*Sx + de*Ex + dr*Rx + di*Iy + ds*Sy + de*Ey + dr*Ry mortality <- ms*S + me*E + mi*I + mr*R + ms*Sx + me*Ex + mi*Ix + mr*Rx + mz*Z + ms*Sy + me*Ey + mi*Iy + mr*Ry #NB this is all mortality (infected and not) H <-discharge + mortality # hospitalisation total (replaces all empty beds with new patients) # find out how many admissions we can take for a single room vs a general admission # Note: this makes recruitment weierd because H_sr doesnt have a multiplier but H_ga does. to_import = min(prop_admitted*P, H) # how many infected want to be admitted? #if(to_import < z_cap + di*Z + n*rsevere*Z + mz*Z - Z){ # if they all fit in the single rooms available then bring them in # H_sr = to_import # to_import = 0 #} #else{ # H_sr = z_cap + di*Z + n*rsevere*Z + mz*Z - Z # else bring as many as possible into single rooms... # to_import = to_import - H_sr #} H_sr =0 H_iga = to_import if(H - H_iga - H_sr > 0){ # if there is space in the hospital bring them all in and then fill rest of space with general admissions H_ga = H - H_iga - H_sr } else{ H_ga = 0 } # print(c(H_ga, H_iga, H_sr, H_ga + H_iga + H_sr, H)) # Function for single-room isolation capacity ctr = contact_testing_rate( i*E,S+E, ct) ctrHCW = contact_testing_rate( i*E,S+E, ctHCW) # Patient Force of infection foi <- beta*((bP2P*I+e*bP2P*E)/Nh + (bH2P*Ihcw+e*bH2P*Ehcw)/Nhcw) # HCW Force of infection foiHCW <- beta*((bP2H*I+e*bP2H*E+z*bP2H*Z)/Nh + c*(bP2H*Ix+e*bP2H*Ex)/Nx + y*(bP2H*Iy+e*bP2H*Ey)/Ny + (bH2H*Ihcw+e*bH2H*Ehcw)/Nhcw) # Suspected patient Force of infection foiX <-(beta*c*((bP2P*Ix+e*bP2P*Ex)/Nx + (bH2P*Ihcw+e*bH2P*Ehcw)/Nhcw)) # Tested -ive cohort Force of infection foiY <-(beta*y*((bP2P*Iy+e*bP2P*Ey)/Ny + (bH2P*Ihcw+e*bH2P*Ehcw)/Nhcw)) # Derivatives dS <- - ms*S - foi*S + (1-atr)*hs*H_ga - ds*S - ctr*S dE <- - me*E + foi*S + (1-atrE)*he*H_ga - de*E - ctr*E - i*E - rhealth*E dI <- - mi*I - rsevere*I - tr*I + i*E - di*I dR <- - mr*R - dr*R + rsevere*I + n*rsevere*Z + rhealth*E + (1-atr)*hr*H_ga - ctr*R dSx <- - ms*Sx - foiX*Sx + ctr*S + ctr*Sy - n*Sx + atr*hs*H_ga - ds*Sx dEx <- - me*Ex + foiX*Sx - i*Ex - de*Ex + ctr*E + ctr*Ey + atrE*he*H_ga - n*Ex dIx <- - mi*Ix + i*Ex - di*Ix + tr*I + H_iga - n*Ix + tr*Iy dRx <- - mr*Rx + ctr*R + ctr*Ry - n*Rx + atr*hr*H_ga - dr*Rx dSy <- - ms*Sy - foiY*Sy + n*Sx - ds*Sy - ctr*Sy dEy <- - me*Ey + foiY*Sy - i*Ey - de*Ey - rhealth*Ey - ctr*Ey dIy <- - mi*Iy + i*Ey - di*Iy - tr*Iy - rsevere*Iy dRy <- - mr*Ry + n*Rx - dr*Ry + rsevere*Iy + rhealth*Ey - ctr*Ry dZ <- - mz*Z - di*Z - n*rsevere*Z + n*Ex + n*Ix dP <- 0 #0.129*((population_prevalence*N) - P)*P/(population_prevalence*N) #-H_ga*(hs + he + hr) - H_sr + ds*S + de*E + di*I + di*Ix + dr*R + ms*S + me*E + mi*I + mr*R + ms*Sx + me*Ex + mi*Ix + mr*Rx dTot <- dS + dE + dI + dR + dSx + dEx + dIx + dRx + dZ + dP + dSy + dEy + dIy + dRy dTot_hosp <- dS + dE + dI + dR + dSx + dEx + dIx + dRx + dZ + dSy + dEy + dIy + dRy dTot_x <- dSx + dEx + dIx + dRx dTot_y <- dSy + dEy + dIy + dRy dShcw <- - foiHCW*Shcw - population_prevalence*Shcw dEhcw <- foiHCW*Shcw - i*Ehcw - ctrHCW*Ehcw - rhealth*Ehcw dIhcw <- i*(1-a)*Ehcw - rhealth*Ihcw - tHCW*Ihcw + population_prevalence*Shcw dRhcw <- rhealth*Ihcw + rabsent*Ahcw + rhealth*Ehcw dAhcw <- i*a*Ehcw + ctrHCW*Ehcw + tHCW*Ihcw - rabsent*Ahcw dTot_hcw <- dShcw + dEhcw + dIhcw + dRhcw Inc <- foi*S + foiX*Sx + foiY*Sy # infected_total <- foi*S + he*H_ga + foiX*Sx + H_sr # infected_identified <- foiX*Sx + ctr*E + tr*I + atr*he*H_ga + H_sr infected_total <- foi*S + foiX*Sx + foiY*Sy infected_identified <- ctr*E + tr*I + tr*Iy # EA:not sure about this- it includes some community-acquired cases too ??? admissions.S = hs*H_ga admissions.E = he*H_ga admissions.I = H_sr + H_iga P2P = beta*S*(bP2P*I+e*bP2P*E)/Nh + Sx*beta*c*((bP2P*Ix+e*bP2P*Ex)/Nx) + Sy*beta*y*((bP2P*Iy+e*bP2P*Ey)/Ny) P2H = beta*Shcw*((bP2H*I+e*bP2H*E+z*bP2H*Z)/Nh + c*(bP2H*Ix+e*bP2H*Ex)/Nx + y*(bP2H*Iy+e*bP2H*Ey)/Ny) H2P = beta*S*(bH2P*Ihcw+e*bH2P*Ehcw)/Nhcw + beta*Sx*c*(bH2P*Ihcw+e*bH2P*Ehcw)/Nhcw + beta*Sy*y*(bH2P*Ihcw+e*bH2P*Ehcw)/Nhcw H2H = beta*Shcw*(bH2H*Ihcw+e*bH2H*Ehcw)/Nhcw # Return values output<-c(c(dS,dE,dI,dR,dSx,dEx,dIx,dRx,dSy,dEy,dIy,dRy,dZ,dP,dTot,dTot_hosp,dShcw,dEhcw,dIhcw,dRhcw,dAhcw,dTot_hcw,dTot_x,dTot_y), Inc/Nh*10000, Inc, infected_total, infected_identified, admissions.S, admissions.E, admissions.I, P2P, P2H, H2P, H2H, foiHCW, population_prevalence, prop_admitted = prop_admitted/hi) #names(ics) <- c(c("dS","dE","dI","dR","dSx","dEx","dIx","dRx","dY","dZ","dP","dTot","dTot_hosp","dShcw","dEhcw","dIhcw","dRhcw","dAhcw","dTot_hcw","dTot_x"), "Inc/Nh*10000", "Inc", "infected_total", "infected_identified", "admissions.S", "admissions.E", "admissions.I", "P2P", "P2H", "H2P", "H2H", "foiHCW", "population_prevalence", "prop_admitted") #list(output) list(c(dS,dE,dI,dR,dSx,dEx,dIx,dRx,dSy,dEy,dIy,dRy,dZ,dP,dTot,dTot_hosp,dShcw,dEhcw,dIhcw,dRhcw,dAhcw,dTot_hcw,dTot_x,dTot_y), Inc/Nh*10000, Inc, infected_total, infected_identified, admissions.S, admissions.E, admissions.I, P2P, P2H, H2P, H2H, foiHCW, population_prevalence, prop_admitted = prop_admitted/hi)}# Function for cumulative infected calculationCumulative <- function(incidence,E.t0,I.t0) { totalInf <-0 for (i in 1:length(incidence)) { if (i==1) { totalInf[i] <- incidence[i] + E.t0 + I.t0 } else { totalInf[i] <- totalInf[i-1] + incidence[i] } } return(totalInf)}# Function for calculating contact testing rate:# infected - number of newly infected patients in this time step# susceptible - number of remaining uninfected patients (S+E in model)# n - number of contacts tested e.g. for cruise ship assume 2 per room so 1 test per infectioncontact_testing_rate <- function(infected, susceptible, n){ x= n*infected/susceptible return(x)}calculate_incidence = function(time, ip) { return(min(ip*(2^(floor(c((time - 5))/2.6))), 3*(10^-4)))}testing_rate_hcw = function(time, period_to_test, tHCW){ if(period_to_test != -1){ if(time %% period_to_test ==0){ tHCW = 1 } else{ tHCW = 0 } } return(tHCW)}absent_rate_hcw = function(time, period_to_test, a){ if(period_to_test != -1){ if(time %% period_to_test ==0){ a = 1 } else{ a= a } } return(a)} ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download