Hα emission in the outskirts of galaxies at z=0.4

Rhythm Shimakawa1^∗NAOJ Fellow
<EMAIL_ADDRESS>(RS)

Masayuki Tanaka1
Satoshi Kikuta2
Masao Hayashi1

1National Astronomical Observatory of Japan (NAOJ), National Institutes of Natural Sciences, Osawa, Mitaka, Tokyo 181-8588, Japan
2Center for Computational Sciences, University of Tsukuba, 1-1-1 Tennodai, Tsukuba, Ibaraki 305-8577, Japan

galaxies: general — galaxies: formation — galaxies: evolution — galaxies: halos

This paper reports detections of H$\alpha$ emission and stellar continuum out to approximately 30 physical kpc, and H$\alpha$ directionality in the outskirts of H$\alpha$-emitting galaxies (H$\alpha$ emitters) at $z=0.4$.
This research adopts narrow-band selected H$\alpha$ emitters at $z=0.4$ from the emission-line object catalog by Hayashi et al., 2020, which is based on data in the Deep and Ultradeep layers of the Hyper Suprime-Cam Subaru Strategic Program.
Deep narrow- and broad-band images of 8625 H$\alpha$ emitters across 16.8 deg$^2$ enable us to construct deep composite emission-line and continuum images.
The stacked images show diffuse H$\alpha$ emission (down to $\sim5\times10^{-20}$ erg s$^{-1}$cm$^{-2}$arcsec$^{-2}$) and stellar continuum (down to $\sim5\times10^{-22}$ erg s$^{-1}$cm$^{-2}$Å$^{-1}$arcsec$^{-2}$), extending beyond 10 kpc at stellar masses $>10^9$ M$_\odot$, parts of which may originate from stellar halos.
Those radial profiles are broadly consistent with each other.
In addition, we obtain a dependence of the H$\alpha$ emission on the position angle because relatively higher H$\alpha$ equivalent width has been detected along the minor-axis towards galaxy disks.
While the H$\alpha$ directionality could be attributed to biconical outflows, further research with hydrodynamic simulations is highly demanded to pin down the exact cause.

§ INTRODUCTION

Galaxy outskirts record mass-assembly histories of merger and accretion events in a hierarchical universe [Searle and Zinn, 1978, White and Rees, 1978], and also feedback processes taking place therein [Dekel and Silk, 1986, Heckman et al., 1990].
Past deep wide-field imaging uncovered various diffuse stellar halos of nearby galaxies and hence diverse buildup histories of outer disks [Courteau et al., 2011, Deason et al., 2014, Merritt et al., 2016, Harmsen et al., 2017], which correlate with their physical properties [Elias et al., 2018].
Combined analyses with simulations in the cosmological framework of Lambda Cold Dark Matter claim that, in general, high fractions of inner and outer stellar halos are formed by in-situ and ex-situ stellar components, respectively [Zhang et al., 2018, Merritt et al., 2020, Font et al., 2020].
Furthermore, the diffuse ionized gas and its physical states in the outskirt of nearby galaxies have been well investigated for understanding feedback mechanisms (e.g., [Veilleux et al., 2005] and references therein).
Zhang et al., 2016, Zhang et al., 2018 successfully detected H$\alpha+$[Nii] emission at $z=$ 0.05–0.2 out to 100 kpc by combining millions of fiber spectra from the Sloan Digital Sky Survey, providing a unique insight into cool gas from associated halos.
However, it remains significantly challenging to observe the galaxy outskirts beyond the local universe due to cosmological surface brightness dimming, which has hampered direct observations to date from the earlier phase of the mass assembling.

Nowadays, a giant wide-field imager on the Subaru Telescope with an 8.2 m primary mirror, Hyper Suprime-Cam [Miyazaki et al., 2018, Furusawa et al., 2018, Kawanomoto et al., 2018, Komiyama et al., 2018] has enabled very deep wide-field observations under seeing conditions of $<1$ arcsec that provides new insights into the outskirts of nearby galaxies [Fukushima et al., 2019, Okamoto et al., 2019].
In addition, such high-quality imaging data over $>100$ deg$^2$ taken by Hyper Suprime-Cam Subaru Strategic Program (HSC-SSP; [Aihara et al., 2018]) can establish even deeper images by stacking numerous objects.
This enables us to investigate the outskirts of massive galaxies beyond the local universe [Huang et al., 2018, Wang et al., 2019], and their environmental dependence [Huang et al., 2018].

With such scientific backgrounds, this research aims at investigating the outskirts of H$\alpha$ line emission from star-forming galaxies at $z=0.4$.
H$\alpha$ line is used extensively as a tracer of star formation [Kennicutt, 1998], which helps us study galaxy outskirts from a unique standpoint.
Based on the emission-line object catalog available from the Second Public Data Release of HSC-SSP (HSC-SSP PDR2; [Aihara et al., 2019, Hayashi et al., 2020]), we perform a stacking analysis of H$\alpha$ line and continuum images for 8625 H$\alpha$ emitting galaxies termed as H$\alpha$ emitters at $z=0.4$.
A wide and deep narrow-band survey by HSC-SSP successfully detects H$\alpha$ emitters at $z=$0.393–0.404 down to $\sim1\times10^{-17}$ erg s$^{-1}$cm$^{-2}$ over 16.8 deg$^2$ (section <ref>).
We investigate radial profiles of their composite H$\alpha$ and continuum images in different stellar mass bins, and study star formation in the outskirts of H$\alpha$ emitters (section <ref>).
We also discuss other possible contributions to H$\alpha$ emission on the galaxy outskirts by shape-aligned stacked images, and lastly summarize the entire flow of this work (section <ref>).

This research adopts the AB magnitude system [Oke and Gunn, 1983] and a Chabrier, 2003 stellar initial mass function (IMF).
Also, we assume cosmological parameters of $\Omega_M=0.310$, $\Omega_\Lambda=0.689$, and $H_0=67.7$ km s$^{-1}$Mpc$^{-1}$ in a flat Lambda cold dark matter model, which are consistent with those from the Planck 2018 VI results [Planck Collaboration et al., 2020].

§ DATA AND METHODOLOGY

The source catalog and data set underlying this paper are available on the HSC-SSP Public Data Release site[<https://hsc.mtk.nao.ac.jp/ssp/data-release/>].
This section briefly overviews the narrow-band emitter catalog in HSC-SSP PDR2 [Hayashi et al., 2020] and describes our procedure of image stacking of broad-band and narrow-band data for the targets, H$\alpha$ emitters at $z=0.4$.

§.§ HSC-SSP PDR2

Narrow-band (H$\alpha$+[Nii]) flux versus logarithmic stellar mass with stellar mass distributions on the top panel.
The light-gray, gray, and black symbols indicate unused emitters in this work, those for the main analysis (section <ref>), and emitters with high ellipticities ($\epsilon>0.33$) that are used in the discussion section (see section <ref>).
The horizontal line depicts the selection threshold for the narrow-band flux ($4\times10^{-17}$ erg s$^{-1}$cm$^{-2}$).
The dotted vertical lines represent selection cuts into three stellar-mass bins.

Summary of the sample selection from the original source catalog ([Hayashi et al., 2020]; $N=8625$).
See also figure <ref> for the selection thresholds.

Sample            Limits      N

All           $f_\mathrm{NB}$$>$4$\times$10$^{-17}$, FWHM$<$0.8     6324

Low-mass      M$_\star$=$10^8$-$10^9$           2121

Mid-mass      M$_\star$=$10^9$-$10^{10}$        3139

High-mass     M$_\star$=$10^{10}$-$10^{11}$      948

Shape-aligned     M$_\star$=$10^9$-$10^{11}$, $\epsilon$$>$0.33     1430

This work is based on the H$\alpha$ emitter sample at $z=0.4$ established by Hayashi et al., 2020; they reported narrow-band selected emission-line objects at $z<2$ based on the Deep and Ultradeep layers in HSC-SSP PDR2 [Aihara et al., 2019].
The catalog contains 8625 H$\alpha$ emitters down to a limiting flux of $\sim1\times10^{-17}$ erg s$^{-1}$cm$^{-2}$ (figure <ref>) and line equivalent width (EW$_\mathrm{H\alpha+[N{\scriptscriptstyle II}]}$) of $>25$ Å in the rest-frame <cit.>.
They show flux excesses in narrow-band, NB921 ($\lambda_\mathrm{center}=9205$ Å), relative to $z$ and $y$-bands by capturing H$\alpha$ ([Nii] doublet also falls into the filter, however, we omit it hereafter) emission line at $z=$0.393–0.404 <cit.>.
The survey area covers 22.09 deg$^2$ (16.79 deg$^2$ excluding bright star mask regions) split into 4 fields (Extended COSMOS, SXDS, ELAIS-N1, DEEP2-3; <cit.>).
Stellar masses of our sample are derived by Mizuki, a SED-based photo-$z$ code [Tanaka, 2015, Tanaka et al., 2018], by fixing the source redshift of $z=0.4$ (stellar_mass_mizuki_zfixed_convflux in the source catalog; [Hayashi et al., 2020]).

To collect H$\alpha$ emitters homogeneously across the entire survey area, we set several selection criteria as follows. We adopt only H$\alpha$ emitters with narrow-band fluxes greater than $4\times10^{-17}$ erg s$^{-1}$cm$^{-2}$ (SFR $=0.1$ M$_\odot$yr$^{-1}$ without dust correction), where the imaging depths achieve approximately 100 percent completeness over the whole survey fields [Hayashi et al., 2020].
Also, this study does not consider areas taken under relatively bad seeing conditions (FWHM $>0.8$ arcsec) in either of narrow-band (NB921) or $z$-band. We then select 6324 H$\alpha$ emitters through these thresholds.

§.§ Stacking analysis

We perform median stacking for narrow-band (NB921) and $z$-band images of the selected H$\alpha$ emitters at $z=0.4$ to derive their typical radial profiles of H$\alpha$ and stellar surface densities.
For the stacking analysis, we obtain coadd images in the NB921 and $z$-band of H$\alpha$ emitters from the Third Public Release of HSC-SSP [Aihara et al., 2021][<https://hsc-release.mtk.nao.ac.jp/das_cutout/pdr3/> (registration is required)].
The data were generated by the dedicated reduction pipeline (hscPipe version 8; [Bosch et al., 2018]). Their spatial resolutions are matched to FWHM $=0.8$ arcsec by Gaussian smoothing.
We confirm that the composite PSF distributions of the NB921 and $z$-band are consistent within a margin of less than 6 percent, which is negligibly small for this study.
Consequently, we derived median pixel values of the target images in each filter.
Here, we separate the emitter sample into three stellar mass bins ($10^{8\mbox{-}9}$, $10^{9\mbox{-}10}$,$10^{10\mbox{-}11}$ M$_\odot$; see figure <ref> and table <ref>).
Based on the composite NB921 and $z$ images, we establish emission line and stellar continuum maps through the following calculation.

\begin{eqnarray}
F_\mathrm{H\alpha} &=& \Delta_\mathrm{NB} \frac{f_\mathrm{NB}-f_\mathrm{BB}}{1-\Delta_\mathrm{NB}/\Delta_\mathrm{BB}},\\
f_\mathrm{c} &=& \frac{f_\mathrm{BB}-f_\mathrm{NB}\cdot\Delta_{NB}/\Delta_\mathrm{BB}}{1-\Delta_\mathrm{NB}/\Delta_\mathrm{BB}},
\end{eqnarray}

where $f_\mathrm{NB}$ and $f_\mathrm{BB}$ are flux densities of the NB921 and $z$-band, and $\Delta_\mathrm{NB}$ and $\Delta_\mathrm{BB}$ indicate filter band widths of the NB921 and $z$-band, respectively [Hayashi et al., 2020].
Furthermore, we adopt a single color-term correction of $f_{BB}=1.04\times f_{z}$ to correct a difference of the central filter wavelength between the NB921 and $z$-band (296 Å), which is consistent with a typical value of high-mass H$\alpha$ emitters.
Although the fixed color-term correction does not significantly affect our main result, this is an important issue on the discussion part regarding to the radial profile of EW$_\mathrm{H\alpha}$ (see section <ref> and Appendix <ref>).

Besides, we generate random 100 realizations of composite emission-line and stellar continuum images to measure flux uncertainties of the stacked images using a bootstrap re-sampling from the same samples but with replacement.
Based on the 100 random combined images, we calculate a standard deviation of self flux counts at a given radius.
We also derive a background deviation given a pixel area ($S$) based on randomly-positioned empty apertures.
Obtained background deviations correlate with the background areas by $\propto S^{\sim0.8}$, suggesting partial pixel-to-pixel correlations (cf. no correlation and full correlation if $\propto S^{0.5}$ and $\propto S$, respectively).
This research defines the mean square errors of deviations of the self-flux and background noises as a total uncertainty of derived flux surface density (section <ref>).
In addition, residual background counts are sampled by measuring median pixel values in the backgrounds of the randomly combined images.
Obtained residual background counts are $\sim4\times10^{-4}$ and $\sim4\times10^{-3}$ in H$\alpha$ and continuum images, respectively, in the zero point magnitude of 27 mag.
These residual backgrounds are subtracted when deriving the flux densities.

§ RESULTS

Radial profiles of surface brightness of H$\alpha$ line emission (blue squares) and normalized stellar continuum (yellow circles) for H$\alpha$ emitters at $z=0.4$.
Less than $2\sigma$ detections are represented by the open symbols.
The error-bars depict $1\sigma$ errors and the light-blue regions are $1\sigma$ errors of the curve fitting (see text).
From the left to right, figures depict radial profiles for H$\alpha$ emitters with stellar masses of (a) $10^{8\mbox{-}9}$, (b) $10^{9\mbox{-}10}$, and (c) $10^{10\mbox{-}11}$ M$_\odot$, respectively.
The dotted curves indicate scaled radial profiles of the composite PSF (fitted by a Moffat distribution).
The dot-dashed curves are the best-fit Gaussian$+$exponential radial profiles, which are expected to trace bulge$+$disk components of the H$\alpha$ emitters.

Median radial profiles of H$\alpha$ emission ($\Sigma F_\mathrm{H\alpha}$) and stellar continuum in three different stellar masses are shown in figure <ref>, where radial continuum distributions are normalized at the peak surface densities in H$\alpha$ line.
Examples of the composite H$\alpha$ and continuum images are represented in Appendix <ref>.
As a result, we detect both H$\alpha$ and continuum of the sample in outer regions of H$\alpha$ emitters at $z=0.4$, especially at the highest stellar mass bin ($10^{10\mbox{-}11}$ M$_\odot$).
While we observe differentials between H$\alpha$ and stellar distributions in $\sim10$ ph-kpc at the highest stellar mass as explored in more detail in section <ref>, they are generally consistent with each other.

To investigate stellar halo contributions, we fit the radial profiles with a combined function of Gaussian ($G$), exponential ($E$), and power law ($P$) distributions by assuming that they respectively trace galactic bulge, disk, and stellar halo:

\begin{eqnarray}
G(r) &=& A_g \exp(-r^2 / w_g),\\
E(r) &=& A_e \exp(-r / r_s),\\
P(r) &=& A_p r^\gamma,
\end{eqnarray}

where $r$ is the physical radius (physical kpc, hereafter ph-kpc) from the center; $A_{g}$, $A_{e}$, $A_{p}$, $w_g$, $r_s$, and $\gamma$ are fitting parameters.
We do not adopt the de Vaucouleurs-law profile [de Vaucouleurs, 1948] for bulge components because our targets are barely resolved or unresolved in the seeing-limited data (FWHM $=4.4$ ph-kpc). We perform a chi-square fitting with the model function shown above using a curve fitting code, lmfit [Newville et al., 2014].

Best-fit parameters and $1\sigma$ errors in three different stellar mass bins
(see figure <ref>).

Mass range      $10^8$–$10^9$     $10^9$–$10^{10}$     $10^{10}$–$10^{11}$

1c($N=$)     (2121)     (3139)     (948)

$A_{g}$ ($\times10^{-17}$)     $3.0\pm0.2$      $4.9\pm0.2$      $14.2\pm0.6$

$w_{g}$                        $11.6\pm0.5$     $15.9\pm0.5$     $20.9\pm0.8$

$A_{e}$ ($\times10^{-17}$)     $1.3\pm0.4$      $1.7\pm0.4$      $3.7\pm1.0$

$r_{s}$                        $3.1\pm0.3$      $3.8\pm0.3$      $4.2\pm0.4$

$A_{p}$ ($\times10^{-18}$)     $0.6\pm2.1$      $1.5\pm2.6$      $11.8\pm9.1$

$\gamma$                       $-1.2\pm1.3$     $-1.2\pm0.5$     $-1.4\pm0.2$

The best-fit parameters and associated errors for the H$\alpha$ flux surface densities are summarized in table <ref>.
Figure <ref> shows that there are residual components beyond $\gtrsim15$ ph-kpc compared to the best-fit bulge$+$disk profiles in the intermediate and high stellar mass bins (figure <ref>b,c).
Such residuals have been commonly observed in stellar distributions of nearby galaxies (e.g., [Courteau et al., 2011, Merritt et al., 2016]), suggesting that they could originate respectively from in-situ hot young stars and old stars in outer disks and stellar halos (see section <ref> for more discussion).
The estimated power law slopes in the outskirts, $\gamma=-1.4\sim-1.2$, are flatter than those of most nearby galaxies $\gamma<-2$ [Harmsen et al., 2017].
However, given the low signal detection, it may be too early to conclude that H$\alpha$ emitters at $z=0.4$ have much shallower power-law slopes than galaxies in the local universe.
Recent studies have reported that stellar halo distributions can be fitted by broken power law; the outer slopes become significantly steeper at broken radii of a few tens of ph-kpc [Deason et al., 2014].
This may be one of the reasons of non detection in both H$\alpha$ and continuum images at $r>30$ ph-kpc.

Moreover, we convert the H$\alpha$ surface densities to SFR surface densities ($\Sigma\mathrm{SFR}$) and investigate fractions of star formation in stellar halos at $>10$ ph-kpc (figure <ref>).
We follow the Kennicutt, 1998 calibration and scale by $1/1.7$ to apply the Chabrier, 2003 IMF.
Also, we assume 20 percent contributions of [Nii]$\lambda\lambda6548,6583$ lines to narrow-band fluxes.
Here it should be noted that we assume no extinction and do not consider warm ionized medium [Reynolds et al., 1973, Haffner et al., 2009].
Thus, the derived $\Sigma\mathrm{SFR}$ values should have significant uncertainties beyond the measurement errors.
In particular, H$\alpha$ extinction in the central components of emitters would exceed 1 mag at the high stellar mass bin [Sobral et al., 2016].
However, central star formation is not the scope of this paper; therefore, we solely focus on the outskirts of H$\alpha$ emitters.
The stacked data have reached $\Sigma\mathrm{SFR}$ down to $\sim3\times10^{-6}$ M$_\odot$yr$^{-1}$kpc$^{-2}$, about a half of which appears to originate from stellar halos.
As inferred from the best-fit results, major fluxes beyond 30 ph-kpc are dominated by halo components; however, significant uncertainties are associated with these results.

(Top panel) Radial profiles of surface SFR densities of H$\alpha$ emitters in three stellar mass bins as in figure <ref> and table <ref> (blue: $10^{8\mbox{-}9}$, green $10^{9\mbox{-}10}$, and red: $10^{10\mbox{-}11}$ M$_\odot$, respectively).
(Bottom panel) Fractions of power law components to the best-fitted radial profiles ($f_\mathrm{halo}$) in respective stellar mass bins at $r>10$ ph-kpc.
Non detection regions with $<2\sigma$ are extrapolated by the best-fitted curves (dashed curves).

§ DISCUSSION AND SUMMARY

This work has identified diffuse H$\alpha$ emission and stellar continuum beyond 10 ph-kpc and up to 30 ph-kpc for the H$\alpha$ emitters at $z=0.4$ based on the deep narrow-band and broad-band stacking taken from HSC-SSP.
There are residual excesses in both H$\alpha$ line and stellar continuum compared to exponential distributions at $\gtrsim20$ ph-kpc.
While such residual emissions could originate respectively from in-situ hot young stars and old stars in outer disks and stellar halos, we cannot ignore contributions from associated halos [Zhang et al., 2016, Zhang et al., 2018] and escape of ionizing radiation (e.g., [Leitherer et al., 1995, Heckman et al., 2011]).
Thus physical origins of diffuse emission lines are still controversial.
In addition, we do not obtain a clear discrepancy between the radial profiles of the emission line and the continuum on the outskirts.
The observed H$\alpha$ surface densities monotonically decline in the outskirts with flatter power-law slopes of $\gamma=-1.4\sim-1.2$ than those in nearby galaxies; however, we note weak detection of some diffuse stellar halos.

Lastly, we delve into causal factors that may contribute to the observed diffuse H$\alpha$ emission on the outskirt by constructing shape-aligned composite images of H$\alpha$ emitters.
This is originally motivated to observe diffuse ionized gas arisen from biconical outflows from galaxy disks as typified by nearby starburst, M82 (e.g., [Lynds and Sandage, 1963, Visvanathan and Sandage, 1972, Bland and Tully, 1988, Veilleux et al., 2003, Yoshida et al., 2019]).
Based on the second moments of the object intensity termed adaptive moments, we calculate ellipticities ($\epsilon$) and position angles of the targets [Bernstein and Jarvis, 2002].
We use the adaptive moments in the $i$-band (i_sdssshape_shape) because the $i$-band data have better and more homogeneous seeing sizes ($\sim0.6$ arcsec) over the survey field.
Also, we employ only elliptic sources at M$_\star>10^{9}$ M$_\odot$ with $\epsilon>0.33$ (i.e., major-axis $>1.5\times$ minor-axis) to assure the credibility of position angle estimates.
We align objects along the major-axis; consequently, we conduct the median stacking following a similar method described in section <ref> for 1430 objects that satisfy these selection criteria (see figure <ref> and table <ref>).

Figure <ref> represents radial profiles of the rest-frame EW$_\mathrm{H\alpha}$ along $\pm30$ deg of the major- and minor-axis.
For deriving radial EW$_\mathrm{H\alpha}$ profiles appropriately, the color term correction is important because the color term is significantly different between the major- and minor-axis.
We thus correct the color term effect in each radial bin on the best effort basis, by using the control sample (see Appendix <ref> for the detailed procedure).
The result shows a clear directivity of EW$_\mathrm{H\alpha}$; for instance, EW$_\mathrm{H\alpha}$ is higher on the outskirts along the minor-axis of the stellar disk (figure <ref>).
A possible factor producing such high EW$_\mathrm{H\alpha}$ could be biconical galactic outflows.
However, as discussed above, there are various potential contributions such as star formation in outer disks and emission from associated halos; therefore, the actual cause would be more complex and complicated.
On the other hand, the moderate dip in the center can be explained by higher nebular dust extinction and/or lower specific star formation rate in the galactic bulge.

In summary, we have unveiled stellar and H$\alpha$ radial profiles on the outskirts of H$\alpha$ emitters at $z=0.4$ for the first time.
However, contributors to diffuse H$\alpha$ emission remain uncertain.
Comparing our results with hydrodynamic simulations and absorption-line observations (e.g., [Werk et al., 2014, Prochaska et al., 2017]) will help us resolve its causal factors from theoretical and independent points of view, respectively.
In addition, more detailed investigations of physical properties of diffuse halos, e.g., derivation of stellar populations on the outskirts based on multi-band image stacking, will deliver further insights into mass-assembly histories surrounding galaxies beyond the local universe.

Radial profiles of rest-frame EW of H$\alpha$ emission in the directions of $\pm30$ deg from major-axis (orange symbols) and minor-axis (purple symbols), respectively.
An image illustration of the areas used in the stacking analysis is denoted at the upper left corner. The error-bars and color-filled regions represent $1\sigma$ errors.
The open symbols indicate $2\sigma$ upper limits.
The dashed curves depict extrapolated lines of the best-fitted models.
The gray region indicates the standard deviations of non-aligned EW$_\mathrm{H\alpha}$ profile based on 100 composite images with random rotation.

Based on data collected at the Subaru Telescope and retrieved from the HSC data archive system, which is operated by Subaru Telescope and Astronomy Data Center at National Astronomical Observatory of Japan. We are honored and grateful for the opportunity of observing the Universe from Maunakea, which has the cultural, historical and natural significance in Hawaii.

The Hyper Suprime-Cam (HSC) collaboration includes the astronomical communities of Japan and Taiwan, and Princeton University. The HSC instrumentation and software were developed by the National Astronomical Observatory of Japan (NAOJ), the Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU), the University of Tokyo, the High Energy Accelerator Research Organization (KEK), the Academia Sinica Institute for Astronomy and Astrophysics in Taiwan (ASIAA), and Princeton University. Funding was contributed by the FIRST program from Japanese Cabinet Office, the Ministry of Education, Culture, Sports, Science and Technology (MEXT), the Japan Society for the Promotion of Science (JSPS), Japan Science and Technology Agency (JST), the Toray Science Foundation, NAOJ, Kavli IPMU, KEK, ASIAA, and Princeton University.
This paper makes use of software developed for the Large Synoptic Survey Telescope. We thank the LSST Project for making their code available as free software at <http://dm.lsst.org>

The Pan-STARRS1 Surveys (PS1) have been made possible through contributions of the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, the Space Telescope Science Institute, the National Aeronautics and Space Administration under Grant No. NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, the National Science Foundation under Grant No. AST-1238877, the University of Maryland, and Eotvos Lorand University (ELTE) and the Los Alamos National Laboratory.

This research was financially supported by JSPS KAKENHI Grant Number JP19K14766.
We would like to thank Editage (<www.editage.com>) for English language editing.
This work made extensive use of the following tools, NumPy
[Harris et al., 2020], Matplotlib [Hunter, 2007], lmfit
[Newville et al., 2014], the Tool for OPerations on Catalogues And Tables,
TOPCAT [Taylor, 2005], a community-developed core Python package for
Astronomy, Astopy [Astropy Collaboration et al., 2013], and Python Data Analysis Library
pandas [McKinney et al., 2010].

§ RADIAL COLOR TERM DEPENDENCE

This section describes how we determine the color term values between NB921 and $z$-band to investigate the H$\alpha$ directionality for our H$\alpha$ emitters at $z=0.4$ in section <ref>.
The shape-aligned composite H$\alpha$ and continuum images of H$\alpha$ emitters are shown in figure <ref>, where we assume a fixed color term value of $\zeta\equiv f_\mathrm{NB}/f_{z}=1.0375$.
The original emitter source catalog [Hayashi et al., 2020] applied the color-term correction by using weighted combined $zy$ fluxes, in particular, to separate [Oii] emitters at $z=1.47$ from distant red galaxies at $z\sim1.3$.
However, we decided not to use $y$-band data in the stacking analyses taking into account its by $\sim1$ mag shallower imaging depth compared to that in the $z$-band.
In addition, we detect significantly brighter PSF wings in the $y$-band than those in the NB921 and $z$-bands, which make the issue even more complicated.

The left and right panels show the composite H$\alpha$ (+[Nii]) and continuum images of shape-aligned H$\alpha$ emitters (see section <ref>).
The white circles depict the radius of 30 ph-kpc.
Here they are stretched to a logarithmic scale to improve visibility.

Instead, we derive typical color term values by using 37163 non-emitter samples with similar stellar populations (M$_\star=10^{9-11}$ M$_\odot$ and SFR $>0.1$ M$_\odot$yr$^{-1}$) at $z_\mathrm{photo}=$ 0.3–0.5 in the same field.
Their stellar masses, SFRs, and photometric redshifts are obtained from Mizuki SED fitting with reduced chi-squares $<5$ [Tanaka, 2015, Tanaka et al., 2018].
The control samples also satisfy the same selection threshold of the shape measurement as for the shape-aligned H$\alpha$ emitter sample ($\epsilon>0.33$).
We then generate 100 realizations of shape-aligned stacked NB921 and $z$ images of randomly selected 10000 non-emitters from the control samples.
We here adjust the sampling rates of non-emitters to match the stellar mass distribution to the H$\alpha$ emitter sample.
Based on random composite images, we measure the median color term values, $\zeta(r)\equiv f_\mathrm{NB}(r)/f_{z}(r)$, in each radial bin along the major- and minor-axis (figure <ref>).

We adopt measured $\zeta(r)$ values to obtain color-term corrected radial profiles of EW$_\mathrm{H\alpha}$ in the major- and minor-axis (figure <ref>).
First of all, figure <ref> indicates bluer color terms along the major-axis than the minor-axis at $\sim10$ ph-kpc, which should be caused by disk components of star-forming galaxies.
Such a low color term leads to a strange dip of EW$_\mathrm{H\alpha}$ in the major-axis if we do not consider this effect (figure <ref>c).
The color terms then increase at larger radii where old stellar populations are more dominant.
In contrast, they tend to decline at $\gtrsim14$ ph-kpc in the minor-axis and $\gtrsim20$ ph-kpc in the major-axis.
Such depressions are thought to be simply due to lower signals in NB921.
However, one should note that the color term uncertainties in the faint end do not affect our conclusions given large EW$_\mathrm{H\alpha}$ errors (figure <ref>).

(Top panel) Radial profiles of surface flux densities in NB921 ($f_\mathrm{NB}$) and $z$-band ($f_\mathrm{z}$).
The orange squares and circles respectively indicate $f_\mathrm{NB}$ and $f_\mathrm{z}$ in the major-axis.
Similarly, the purple squares and circles show $f_\mathrm{NB}$ and $f_\mathrm{z}$ in the minor-axis.
The black cross and plus symbols depict composite PSFs in NB921 and $z$-band, respectively.
The dotted line is the best-fit Moffat model to the composite PSF in NB921, adopted throughout the paper.
(Bottom panel) Radial color term distribution along the major-axis (orange) and minor-axis (purple).
In both panels, error-bars indicate 68 percentiles in 100 random realizations.

(a) Same as in figure <ref> but for the major-axis of the shape-aligned composite images.
Those with color-term corrections based on $\zeta(r)$ derived in figure <ref>, are shown by blue and yellow open symbols, respectively.
(b) Same as in figure <ref>a, but for the minor-axis.
(c) Same as in figure <ref>, but we here compare radial profiles of EW$_\mathrm{H\alpha}$ with and without color term corrections, which are respectively depicted by open symbols with dark-colored regions and filled symbols with light-colored regions.

[Aihara et al., 2018]
Aihara, H., et al. 2018, , 70, S4
[Aihara et al., 2019]
Aihara, H., et al. 2019, , 71, 114
[Aihara et al., 2021]
Aihara, H., et al. 2021, arXiv e-prints, arXiv:2108.13045
[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, et al. 2013, , 558, A33
[Bernstein and Jarvis, 2002]
Bernstein, G. M., & Jarvis, M. 2002, , 123, 583
[Bland and Tully, 1988]
Bland, J., & Tully, B. 1988, , 334, 43
[Bosch et al., 2018]
Bosch, J., et al. 2018, , 70, S5
[Chabrier, 2003]
Chabrier, G. 2003, , 115, 763
[Courteau et al., 2011]
Courteau, S., et al. 2011, , 739, 20
[Deason et al., 2014]
Deason, A. J., et al. 2014, , 787, 30
[de Vaucouleurs, 1948]
de Vaucouleurs, G. 1948, Annales d'Astrophysique, 11, 247
[Dekel and Silk, 1986]
Dekel, A., & Silk, J. 1986, , 303, 39
[Elias et al., 2018]
Elias, L. M., et al. 2018, , 479, 4004
[Font et al., 2020]
Font, A. S., et al. 2020, , 498, 1765
[Fukushima et al., 2019]
Fukushima, T., et al. 2019, , 71, 72
[Furusawa et al., 2018]
Furusawa, H., et al. 2018, , 70, S3
[Haffner et al., 2009]
Haffner, L. M., et al. 2009, Reviews of Modern Physics, 81, 969
[Harmsen et al., 2017]
Harmsen, B., et al. 2017, , 466, 1491
[Harris et al., 2020]
Harris, C. R., et al. 2020, , 585, 357
[Hayashi et al., 2018]
Hayashi, M., et al. 2018, , 70, S17
[Hayashi et al., 2020]
Hayashi, M., et al. 2020, , 72, 86
[Heckman et al., 1990]
Heckman, T. M., Armus, L., & Miley, G. K. 1990, , 74, 833
[Heckman et al., 2011]
Heckman, T. M., et al. 2011, , 730, 5
[Huang et al., 2018]
Huang, S., et al. 2018, , 475, 3348
[Huang et al., 2018]
Huang, S., et al. 2018, , 480, 521
[Hunter, 2007]
Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90
[Kawanomoto et al., 2018]
Kawanomoto, S., et al. 2018, , 70, 66
[Kennicutt, 1998]
Kennicutt, R. C. 1998, , 36, 189
[Komiyama et al., 2018]
Komiyama, Y., et al. 2018, , 70, S2
[Leitherer et al., 1995]
Leitherer, C., et al. 1995, , 454, L19
[Lynds and Sandage, 1963]
Lynds, C. R., & Sandage, A. R. 1963, , 137, 1005
[McKinney et al., 2010]
McKinney W., et al. 2010, Proc. 9th Python Sci. Conf., 1697900, 51
[Merritt et al., 2016]
Merritt, A., et al. 2016, , 830, 62
[Merritt et al., 2020]
Merritt, A., et al. 2020, , 495, 4570
[Miyazaki et al., 2018]
Miyazaki, S., et al. 2018, , 70, S1
[Newville et al., 2014]
Newville, M., et al. 2014, Zenodo,
[Niederhofer et al., 2016]
Niederhofer, F., et al. 2016, , 592, A47
[Okamoto et al., 2019]
Okamoto, S., et al. 2019, , 884, 128
[Oke and Gunn, 1983]
Oke, J. B., & Gunn, J. E. 1983, , 266, 713
[Reynolds et al., 1973]
Reynolds, R. J., Scherb, F., & Roesler, F. L. 1973, , 185, 869
[Planck Collaboration et al., 2020]
Planck Collaboration, Aghanim, et al. 2020, , 641, A6
[Searle and Zinn, 1978]
Searle, L., & Zinn, R. 1978, , 225, 357
[Sobral et al., 2016]
Sobral, D., et al. 2016, , 458, 3443
[Prochaska et al., 2017]
Prochaska, J. X., et al. 2017, , 837, 169
[Tanaka, 2015]
Tanaka, M. 2015, , 801, 20
[Tanaka et al., 2018]
Tanaka, M., et al. 2018, , 70, S9
[Taylor, 2005]
Taylor, M. B. 2005, Astronomical Data Analysis Software and Systems XIV, 347, 29
[Veilleux et al., 2003]
Veilleux, S., et al. 2003, , 126, 2185
[Veilleux et al., 2005]
Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, , 43, 769
[Visvanathan and Sandage, 1972]
Visvanathan, N., & Sandage, A. 1972, , 176, 57
[Wang et al., 2019]
Wang, W., et al. 2019, , 487, 1580
[Werk et al., 2014]
Werk, J. K., et al. 2014, , 792, 8
[White and Rees, 1978]
White, S. D. M., & Rees, M. J. 1978, , 183, 341
[Yoshida et al., 2019]
Yoshida, M., et al. 2019, , 71, 87
[Zhang et al., 2016]
Zhang, H., et al. 2016, , 833, 276
[Zhang et al., 2018]
Zhang, H., Zaritsky, D., & Behroozi, P. 2018, , 861, 34
[Zhang et al., 2018]
Zhang, J., et al. 2018, , 855, 78

=== Frequently used abbreviation of journal names ===
Phys. Rev. D
Phys. Rev. Lett.
Sol. Phys.
Space Sci. Rev.
IAU Circ.
Geochim. Cosmochim. Acta
J. Geophys. Res.
Mem. Soc. Astron. Italiana
Nucl. Phys. A
Proc. SPIE
AIP Conf. Proc.
ASP Conf. Ser.


# Supporting Information for “Quantifying the HIV reservoir with dilution
assays and deep viral sequencing”

by Sarah C. Lotspeich, Brian D. Richardson, Pedro L. Baldoni,
Kimberly P. Enders and Michael G. Hudgens

##  Web Appendix A Finding the MLE

Note that, for the model assuming perfect assay sensitivity and specificity,
the maximum likelihood estimator (MLE) must satisfy
$\widehat{\lambda}_{i}\in(0,\infty)$ as long as neither of the extreme assay
scenarios from Section 2.3 are met. In this case, an unconstrained
optimization method is used to maximize the log-likelihood with respect to
$\theta_{i}=\log(\lambda_{i})\in\mathbb{R}$, and the resulting maximizers
$\widehat{\theta}_{i}$ are transformed to obtain the MLEs
$\widehat{\lambda}_{i}=\exp(\widehat{\theta}_{i})\in(0,\infty)$. In
particular, the Broyden, Fletcher, Goldfarb, and Shannon (BFGS) procedure, a
quasi-Newton method (bonnans_numerical_2006), is used to find
$\boldsymbol{\widehat{\theta}}=(\widehat{\theta}_{1},\dots,\widehat{\theta_{n}})^{T}$.
This method is implemented via the R (R) function optim using the option
method="BFGS". The initial values are set to
$\boldsymbol{\widehat{\theta}}^{(0)}=\log\\{-\log\left(1-\boldsymbol{Y}/M\right)\\}$,
and the number of iterations is constrained to a maximum of
$\mathrm{10{,}000}$. The gradient of the log-likelihood has the following
closed form,

$\displaystyle\nabla l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})$
$\displaystyle=\begin{bmatrix}\frac{\partial}{\partial\lambda_{1}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})&\dots&\frac{\partial}{\partial\lambda_{n}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\end{bmatrix}^{\textrm{T}}$

with elements given by

$\displaystyle\frac{\partial}{\partial\lambda_{i}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})$
$\displaystyle=\frac{Y_{i}}{\exp\left(\lambda_{i}\right)-1}-\left(M_{N}+m-Y_{i}\right)+\frac{M-(M_{N}+m)}{\exp\left(\Lambda\right)-1},$
(S.1)

for DVL-specific rate $\lambda_{i}$, $i\in\\{1,\dots,n\\}$. Supplying this
gradient to the BFGS routine allows for efficient numerical optimization.

The BFGS procedure is guaranteed to converge to a global maximum if the
objective function is concave (bonnans_numerical_2006). To see that the log-
likelihood is in fact concave, note that its second derivatives are

$\displaystyle\frac{\partial^{2}}{\partial\lambda_{i}^{2}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})=$
$\displaystyle-\frac{Y_{i}\exp\left(\lambda_{i}\right)}{\left\\{\exp\left(\lambda_{i}\right)-1\right\\}^{2}}-\frac{\left\\{M-\left(M_{N}+m\right)\right\\}\exp\left(\Lambda\right)}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{2}},$
(S.2)
$\displaystyle\frac{\partial^{2}}{\partial\lambda_{i}\partial\lambda_{j}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})=$
$\displaystyle-\frac{\left\\{M-\left(M_{N}+m\right)\right\\}\exp\left(\Lambda\right)}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{2}},$
(S.3)

for $i,j\in\\{1,\dots,n\\}$. These derivatives can also be written as

$\displaystyle\frac{\partial^{2}}{\partial\boldsymbol{\lambda}^{2}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})=-$
$\displaystyle\text{diag}\left[\frac{Y_{i}\exp\left(\boldsymbol{\lambda}\right)}{\left\\{\exp\left(\boldsymbol{\lambda}\right)-1\right\\}^{2}}\right]-\frac{\left\\{M-\left(M_{N}+m\right)\right\\}\exp\left(\Lambda\right)}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{2}}\boldsymbol{J}_{n},$

where $\boldsymbol{J}_{n}$ is an $n\times n$ matrix of ones and
$\text{diag}(\boldsymbol{v})$ is an $n\times n$ diagonal matrix with the
$n$-vector $\boldsymbol{v}$ on the diagonal. Since
$\partial^{2}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})/\partial\boldsymbol{\lambda}^{2}$
is a sum of a negative definite matrix and a negative semi-definite matrix, it
is negative definite. Thus the log-likelihood is concave, guaranteeing global
convergence of the BFGS procedure.

As noted in Section 2.3, analytical solutions for
$\widehat{\boldsymbol{\lambda}}$ do exist in two special cases with data from
a single dilution level: (i) when no positive wells are deep sequenced (Web
Appendix A.1) or (ii) when all positive wells are deep sequenced (Web Appendix
A.2).

###  Web Appendix A.1 Closed-form MLE when no deep sequencing information is
available

In case (i), the MLE is equivalent to the IUPM estimator from Myers1994 and
Trumble2017, since without the UDSA data (1) in the main text reduces to the
likelihood from Myers et al. with $D=1$ dilution level, i.e.,

$\displaystyle L(\Lambda|M_{N})$
$\displaystyle=\left\\{1-\exp\left(-\Lambda\right)\right\\}^{(M-M_{N})}\exp\left(-\Lambda\right)^{M_{N}}.$
(S.4)

In fact, in this case the MLE for the IUPM can be found directly, rather than
finding $\widehat{\boldsymbol{\lambda}}$ first and taking its sum. Taking the
derivative of the log-likelihood based on (S.4) with respect to $\Lambda$,

$\displaystyle\frac{\partial}{\partial\Lambda}l(\Lambda|M_{N})$
$\displaystyle=\frac{M-M_{N}}{\exp\left(\Lambda\right)-1}-M_{N},$

and finding its root leads to the MLE:
$\widehat{\Lambda}=\log\left(M/M_{N}\right)$.

###  Web Appendix A.2 Closed-form MLE when deep sequencing information is
fully available

In case (ii), $m=M_{P}$ and therefore (1) has the simplified form

$\displaystyle L(\boldsymbol{\lambda}|\boldsymbol{Y})$
$\displaystyle=\prod_{i=1}^{n}\left\\{1-\exp(-\lambda_{i})\right\\}^{Y_{i}}\exp(-\lambda_{i})^{(M-Y_{i})},$
(S.5)

with corresponding log-likelihood proportional to

$\displaystyle l(\boldsymbol{\lambda}|\boldsymbol{Y})$
$\displaystyle=\sum_{i=1}^{n}Y_{i}\log\left\\{1-\exp\left(-\lambda_{i}\right)\right\\}-\sum_{i=1}^{n}(M-Y_{i})\lambda_{i}.$
(S.6)

The MLE is defined to maximize (S.6); in other words, $\widehat{\lambda}_{i}$
are the roots of the following derivatives with respect to each DVL-specific
rate $\lambda_{i}$:

$\displaystyle\frac{\partial}{\partial\lambda_{i}}l(\boldsymbol{\lambda}|\boldsymbol{Y})$
$\displaystyle=Y_{i}\left\\{\frac{\exp\left(-\lambda_{i}\right)}{1-\exp\left(-\lambda_{i}\right)}\right\\}-(M-Y_{i})$
(S.7)

for $i\in\\{1,\dots,n\\}$. Setting (S.7) equal to zero and solving for
$\lambda_{i}$ yields the MLE:
$\widehat{\boldsymbol{\lambda}}=-\log\left(1-\boldsymbol{Y}/M\right)$.

##  Web Appendix B Derivations for the Variance Estimator

In Section 2.3, it was noted that the $n\times n$ covariance matrix
$\boldsymbol{\Sigma}$ of the MLE $\widehat{\boldsymbol{\lambda}}$ is the
inverse of the Fisher Information, which can be estimated with
$\widehat{\boldsymbol{\Sigma}}=\mathcal{I}(\boldsymbol{\lambda})^{-1}|_{\boldsymbol{\lambda}=\widehat{\boldsymbol{\lambda}}}$.
In this section, $\mathcal{I}(\boldsymbol{\lambda})$ is derived using the log-
likelihood based on (1) from the main text. To begin, consider the case where
data are available from a single dilution level of one million cells per well,
as in Sections 2.1–2.5. In Web Appendix D, an extension to the more general
setting with a different single dilution or multiple dilutions is provided.

To compute the Fisher information, expectations of the random quantities
$Y_{i}$, $M_{N}$, and $m$ are needed. The number of negative wells $M_{N}$ is
a binomial random variable with $M$ trials (i.e., the total number of
replicate wells) and probability of success (i.e., of a well being uninfected)
equal to $\Pr(W_{j}=0)=\exp\left(-\Lambda\right)$, so
$\textrm{E}\left(M_{N}\right)=M\exp\left(-\Lambda\right)$. Given that $M$ and
$q$ are fixed, $m=\lfloor q(M-M_{N})\rceil$ is a function of $M_{N}$ and thus

$\displaystyle\textrm{E}(m)$ $\displaystyle=\sum_{k=0}^{M}\lfloor
q(M-k)\rceil\Pr(M_{N}=k)$ $\displaystyle=\sum_{k=0}^{M}\lfloor
q(M-k)\rceil{M\choose
k}\left\\{\exp(-\Lambda)\right\\}^{k}\left\\{1-\exp(-\Lambda)\right\\}^{M-k}.$

A similar approach can be used to calculate the expectation of $Y_{i}$. Recall
that for DVL $i$, $i\in\\{1,\dots,n\\}$, $Y_{i}=\sum_{j=1}^{M}X_{ij}R_{j}$ is
defined as the total number of wells that are positive for DVL $i$ (i.e.,
$X_{ij}=1$) and have complete data (i.e., $R_{j}=1$). Using this definition
and the assumption that the DVL indicators $X_{ij}$ in the $M$ wells are
independent and identically distributed, the expectation of $Y_{i}$ is
$\textrm{E}(Y_{i})=\sum_{j=1}^{M}\textrm{E}(X_{ij}R_{j})=M\textrm{E}(X_{i1}R_{1}),$
where the first well was chosen here without loss of generality.

Since $X_{i1}R_{1}$ is a product of indicator variables,
$\textrm{E}(X_{i1}R_{1})=\Pr(X_{i1}R_{1}=1)=\Pr(X_{i1}=1,R_{1}=1)$. To compute
this probability, it is useful to condition on $M_{P}$, the number of positive
wells, and employ the Law of Total Probability as follows

$\displaystyle\Pr(X_{i1}=1,R_{1}=1)=\sum_{k=0}^{M}\Pr(X_{i1}=1,R_{1}=1|M_{P}=k)\Pr(M_{P}=k).$
(S.8)

Under the assumption of perfect sensitivity and specificity for the QVOA (A3),
if $X_{i1}=1$ then $W_{1}=1$, also. Therefore,
$\Pr(X_{i1}=1,R_{1}=1|M_{P}=k)=\Pr(W_{1}=1,X_{i1}=1,R_{1}=1|M_{P}=k)$, which
can be factored as follows

$\displaystyle\Pr(W_{1}=1,X_{i1}=1,R_{1}=1|M_{P}=k)$
$\displaystyle=\Pr(R_{1}=1|X_{i1}=1,W_{1}=1,M_{P}=k)\Pr(X_{i1}=1|W_{1}=1,M_{P}=k)\Pr(W_{1}=1|M_{P}=k)$
$\displaystyle=\left(\frac{\lfloor
qk\rceil}{k}\right)\left\\{\frac{1-\exp(-\lambda_{i})}{1-\exp(-\Lambda)}\right\\}\left(\frac{k}{M}\right)$
$\displaystyle=\left(\frac{\lfloor
qk\rceil}{M}\right)\left\\{\frac{1-\exp(-\lambda_{i})}{1-\exp(-\Lambda)}\right\\}.$
(S.9)

Combining (S.8) and (S.9) leads to the expectation of $Y_{i}$:

$\displaystyle\textrm{E}(Y_{i})$ $\displaystyle=M\textrm{E}(X_{i1}R_{1})$
$\displaystyle=M\sum_{k=0}^{M}\Pr(M_{P}=k)\Pr(X_{i1}=1,R_{1}=1|M_{P}=k)$
$\displaystyle=M\sum_{k=1}^{M}{M\choose
k}\left\\{1-\exp(-\Lambda)\right\\}^{k-1}\exp(-\Lambda)^{M-k}\left(\frac{\lfloor
qk\rceil}{M}\right)\left\\{\frac{1-\exp(-\lambda_{i})}{1-\exp(-\Lambda)}\right\\}$
$\displaystyle=\left\\{1-\exp(-\lambda_{i})\right\\}\sum_{k=1}^{M}{M\choose
k}\lfloor
qk\rceil\left\\{1-\exp(-\Lambda)\right\\}^{k-1}\exp(-\Lambda)^{M-k}.$ (S.10)

Since the second derivatives of the log-likelihood are linear in $Y_{i}$,
$M_{N}$, and $m$, the elements of the Fisher information can by computed by
substituting $\textrm{E}(Y_{i})$, $\textrm{E}(M_{N})$, and $\textrm{E}(m)$
into the negations of (S.2) and (S.3) to obtain

$\displaystyle\mathcal{I}_{ii}(\boldsymbol{\lambda})$
$\displaystyle=\textrm{E}\left\\{-\frac{\partial^{2}}{\partial^{2}\lambda_{i}^{2}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\right\\}=\frac{\textrm{E}(Y_{i})\exp\left(\lambda_{i}\right)}{\left\\{\exp\left(\lambda_{i}\right)-1\right\\}^{2}}+\frac{\left\\{M-M\exp(-\Lambda)-\textrm{E}(m)\right\\}\exp\left(\Lambda\right)}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{2}}$
(S.11)

along the diagonal and

$\displaystyle\mathcal{I}_{ij}(\boldsymbol{\lambda})$
$\displaystyle=\textrm{E}\left\\{-\frac{\partial^{2}}{\partial\lambda_{i}\partial\lambda_{j}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\right\\}=\frac{\left\\{M-M\exp(-\Lambda)-\textrm{E}(m)\right\\}\exp\left(\Lambda\right)}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{2}},$
(S.12)

on the off-diagonal, respectively, of $\mathcal{I}(\boldsymbol{\lambda})$.

##  Web Appendix C Derivations for the Bias Correction

In Section 2.4 of the main text, the bias correction for the MLE of the DVL-
specific rates was introduced as
$B(\widehat{\boldsymbol{\lambda}})=\widehat{\boldsymbol{\Sigma}}~{}\boldsymbol{A}(\widehat{\boldsymbol{\lambda}})\mathrm{vec}(\widehat{\boldsymbol{\Sigma}})$,
where
$\boldsymbol{A}\\!\left(\boldsymbol{\lambda}\right)=\left[\boldsymbol{A}_{1}\\!\left(\boldsymbol{\lambda}\right),\dots,\boldsymbol{A}_{n}\\!\left(\boldsymbol{\lambda}\right)\right]$
is the $n\times n^{2}$ matrix with $n\times n$ submatrices

$\displaystyle\boldsymbol{A}_{i}\\!\left(\boldsymbol{\lambda}\right)=\frac{\partial}{\partial\lambda_{i}}\mathcal{I}(\boldsymbol{\lambda})-\frac{1}{2}\textrm{E}\left\\{\frac{\partial^{3}}{\partial\boldsymbol{\lambda}\partial\boldsymbol{\lambda}^{T}\partial\lambda_{i}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\right\\}$

and $\mathrm{vec}\\!\left(\boldsymbol{\Sigma}\right)$ denotes the $n^{2}\times
1$ column vector obtained by stacking the columns of $\boldsymbol{\Sigma}$. In
this section, the necessary components of
$\boldsymbol{A}_{i}(\boldsymbol{\lambda})$ are derived. As with the derivation
of the Fisher information in Web Appendix B, it will be useful to first
consider the case where data are available from one dilution level of one
million cells per well. An extension to the more general setting is provided
in Web Appendix D.

Closed form expressions for the Fisher information are given by (S.11) and
(S.12), and can be differentiated numerically. Since the expectation of
$\partial^{3}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})/\partial\lambda_{i}^{3}$
is used in $\boldsymbol{A}_{i}(\boldsymbol{\lambda})$, an analytic expression
for the third derivative of the log-likelihood is required. Picking up with
the second derivative in (S.2) and (S.3), the third derivative has entries

$\displaystyle\frac{\partial^{3}}{\partial\lambda_{i}^{3}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})$
$\displaystyle=\frac{Y_{i}\exp\left(\lambda_{i}\right)\left\\{\exp\left(\lambda_{i}\right)+1\right\\}}{\left\\{\exp\left(\lambda_{i}\right)-1\right\\}^{3}}+\frac{\left(M-M_{N}-m\right)\exp\left(\Lambda\right)\left\\{\exp\left(\Lambda\right)+1\right\\}}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{3}}$
$\displaystyle\frac{\partial^{3}}{\partial\lambda_{i}^{2}\partial\lambda_{j}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})$
$\displaystyle=\frac{\partial^{3}}{\partial\lambda_{i}\partial\lambda_{j}\partial\lambda_{k}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})=\frac{\left(M-M_{N}-m\right)\exp\left(\Lambda\right)\left\\{\exp\left(\Lambda\right)+1\right\\}}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{3}},$

for distinct $i,j,k\in\\{1,\dots,n\\}$. Then the expectations of the third
derivative are

$\displaystyle\textrm{E}\left\\{\frac{\partial^{3}}{\partial\lambda_{i}^{3}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\right\\}=$
$\displaystyle\frac{\textrm{E}(Y_{i})\exp\left(\lambda_{i}\right)\left\\{\exp\left(\lambda_{i}\right)+1\right\\}}{\left\\{\exp\left(\lambda_{i}\right)-1\right\\}^{3}}$
$\displaystyle+\frac{\left\\{M-M\exp(-\Lambda)-\textrm{E}(m)\right\\}\exp\left(\Lambda\right)\left\\{\exp\left(\Lambda\right)+1\right\\}}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{3}},$
(S.13)
$\displaystyle\textrm{E}\left\\{\frac{\partial^{3}}{\partial\lambda_{i}^{2}\partial\lambda_{j}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\right\\}=$
$\displaystyle\textrm{E}\left\\{\frac{\partial^{3}}{\partial\lambda_{i}\partial\lambda_{j}\partial\lambda_{k}}l(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})\right\\}$
$\displaystyle=$
$\displaystyle\frac{\left\\{M-M\exp(-\Lambda)-\textrm{E}(m)\right\\}\exp\left(\Lambda\right)\left\\{\exp\left(\Lambda\right)+1\right\\}}{\left\\{\exp\left(\Lambda\right)-1\right\\}^{3}},$
(S.14)

where $\textrm{E}(m)$ and $\textrm{E}(Y_{i})$ are derived in Web Appendix B.

##  Web Appendix D Additional Derivations for Multiple Dilution Levels

The derivations in Web Appendices B and C can be extended to the multiple
dilution level setting as described below. Following the notation introduced
in the main text, suppose that assay data
$(M_{N}^{(d)},\boldsymbol{Y}^{(d)},u^{(d)})$, $d\in\\{1,\dots D\\}$, are
available from $D$ distinct dilution levels, $D\in\\{1,2,\dots\\}$. Based on
(4) in the main text, the joint log-likelihood given data from all $D$
dilution levels is proportional to

$\displaystyle\widetilde{l}(\boldsymbol{\tau}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})=\sum_{d=1}^{D}l(u^{(d)}\boldsymbol{\tau}|M_{N}^{(d)},\boldsymbol{Y}^{(d)}),$

where $l(u\boldsymbol{\tau}|M_{N},\boldsymbol{Y})$, based on (1) in the main
text, is proportional to the log-likelihood for a single dilution level.

Using the chain rule, the first two derivatives are

$\displaystyle\frac{\partial}{\partial\boldsymbol{\tau}}\widetilde{l}(\boldsymbol{\tau}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})$
$\displaystyle=\sum_{d=1}^{D}u^{(d)}\frac{\partial}{\partial\boldsymbol{\lambda}}l(u^{(d)}\boldsymbol{\tau}|M_{N}^{(d)},\boldsymbol{Y}^{(d)})$
$\displaystyle\frac{\partial^{2}}{\partial\boldsymbol{\tau}\partial\boldsymbol{\tau}^{\textrm{T}}}\widetilde{l}(\boldsymbol{\tau}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})$
$\displaystyle=\sum_{d=1}^{D}\left(u^{(d)}\right)^{2}\frac{\partial^{2}}{\partial\boldsymbol{\lambda}\partial\boldsymbol{\lambda}^{\textrm{T}}}l(u^{(d)}\boldsymbol{\tau}|M_{N}^{(d)},\boldsymbol{Y}^{(d)}).$

Using the linearity of expectations, the Fisher information matrix
$\widetilde{\mathcal{I}}(\boldsymbol{\tau})$ given data from multiple
dilutions can be expressed in terms of $\mathcal{I}(\boldsymbol{\lambda})$
from Web Appendix B:

$\displaystyle\widetilde{\mathcal{I}}(\boldsymbol{\tau})=\textrm{E}\left\\{-\frac{\partial^{2}}{\partial\boldsymbol{\tau}\partial\boldsymbol{\tau}^{\textrm{T}}}\widetilde{l}(\boldsymbol{\tau}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})\right\\}=\sum_{d=1}^{D}\left(u^{(d)}\right)^{2}\mathcal{I}(u^{(d)}\boldsymbol{\tau}|M_{N}^{(d)},\boldsymbol{Y}^{(d)}),$

where $\mathcal{I}(u^{(d)}\boldsymbol{\tau}|M_{N}^{(d)},\boldsymbol{Y}^{(d)})$
denotes the previously developed Fisher information matrix for
$\boldsymbol{\lambda}$ evaluated at
$(\boldsymbol{\lambda}=u^{(d)}\boldsymbol{\tau},M_{N}=M_{N}^{(d)},\boldsymbol{Y}=\boldsymbol{Y}^{(d)})$.
Then the covariance matrix $\widetilde{\boldsymbol{\Sigma}}$ of
$\boldsymbol{\tau}$ can be estimated with
$\widetilde{\mathcal{I}}^{-1}(\boldsymbol{\tau})|_{\boldsymbol{\tau}=\widehat{\boldsymbol{\tau}}}$.

A bias correction term for the MLE $\widehat{\boldsymbol{\tau}}$ given data
from multiple dilution levels is

$\displaystyle\widetilde{B}\left(\widehat{\boldsymbol{\tau}}\right)=\left\\{\widetilde{\mathcal{I}}^{-1}(\boldsymbol{\tau})|_{\boldsymbol{\tau}=\widehat{\boldsymbol{\tau}}}\right\\}\widetilde{\boldsymbol{A}}\left(\widehat{\boldsymbol{\tau}}\right)\mathrm{vec}\left\\{\widetilde{\mathcal{I}}^{-1}(\boldsymbol{\tau})|_{\boldsymbol{\tau}=\widehat{\boldsymbol{\tau}}}\right\\},$
(S.15)

where
$\widetilde{\boldsymbol{A}}(\boldsymbol{\tau})=[\widetilde{\boldsymbol{A}}_{1}(\boldsymbol{\tau}),\dots,\widetilde{\boldsymbol{A}}_{n}(\boldsymbol{\tau})]$
is the $n\times n^{2}$ matrix with $n\times n$ submatrices defined as

$\displaystyle\widetilde{\boldsymbol{A}}_{i}\left(\boldsymbol{\tau}\right)=\frac{\partial}{\partial\tau_{i}}\widetilde{\mathcal{I}}(\boldsymbol{\tau})-\frac{1}{2}\textrm{E}\left\\{\frac{\partial^{3}}{\partial\boldsymbol{\tau}\partial\boldsymbol{\tau}^{\textrm{T}}\partial\tau_{i}}\widetilde{l}(\boldsymbol{\tau}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})\right\\}.$

As in the single dilution case, the derivative
$\partial\widetilde{\mathcal{I}}(\boldsymbol{\tau})/\partial\tau_{i}$ can be
computed numerically. Also, using the linearity of expectations,

$\displaystyle{\rm
E}\left\\{\frac{\partial^{3}}{\partial\boldsymbol{\tau}\partial\boldsymbol{\tau}^{T}\partial\tau_{i}}\tilde{l}(\widehat{\boldsymbol{\tau}}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})\right\\}=\sum_{d=1}^{D}\left(u^{(d)}\right)^{3}{\rm
E}\left\\{\frac{\partial^{3}}{\partial\boldsymbol{\lambda}\partial\boldsymbol{\lambda}^{T}\partial\lambda_{i}}l(u^{(d)}\widehat{\boldsymbol{\tau}}|M_{N}^{(d)},\boldsymbol{Y}^{(d)})\right\\},$

where the expectations inside the sum on the right-hand side can be computed
using (S.13) and (S.14) from Web Appendix C.

Table S1: Simulation results with a single dilution level and a non-constant
rate of infected cells for all distinct viral lineages. The true overall IUPM
in all settings was $T=1$.

| Without UDSA | With UDSA
---|---|---
| MLE | Bias-Corrected MLE | MLE | Bias-Corrected MLE
$\boldsymbol{n^{\prime}}$ | $\boldsymbol{M}$ | $\boldsymbol{q}$ | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP
6 | 12 | 0.50 | $0.10$ | $0.42$ | $0.45$ | $0.95$ | $0.00$ | $0.42$ | $0.37$ | $0.90$ | $0.09$ | $0.36$ | $0.38$ | $0.94$ | $-0.02$ | $0.36$ | $0.34$ | $0.98$
|  | 0.75 | $0.10$ | $0.42$ | $0.45$ | $0.95$ | $0.00$ | $0.42$ | $0.37$ | $0.90$ | $0.07$ | $0.34$ | $0.35$ | $0.94$ | $0.00$ | $0.34$ | $0.33$ | $0.97$
|  | 1.00 | $0.10$ | $0.42$ | $0.45$ | $0.95$ | $0.00$ | $0.42$ | $0.37$ | $0.90$ | $0.05$ | $0.32$ | $0.33$ | $0.95$ | $0.00$ | $0.32$ | $0.31$ | $0.97$
| 24 | 0.50 | $0.05$ | $0.28$ | $0.30$ | $0.96$ | $0.01$ | $0.28$ | $0.28$ | $0.94$ | $0.04$ | $0.25$ | $0.25$ | $0.95$ | $-0.01$ | $0.25$ | $0.24$ | $0.97$
|  | 0.75 | $0.05$ | $0.28$ | $0.30$ | $0.96$ | $0.01$ | $0.28$ | $0.28$ | $0.94$ | $0.03$ | $0.23$ | $0.24$ | $0.94$ | $0.00$ | $0.23$ | $0.23$ | $0.95$
|  | 1.00 | $0.05$ | $0.28$ | $0.30$ | $0.96$ | $0.01$ | $0.28$ | $0.28$ | $0.94$ | $0.02$ | $0.22$ | $0.23$ | $0.95$ | $0.00$ | $0.22$ | $0.22$ | $0.96$
| 32 | 0.50 | $0.03$ | $0.24$ | $0.25$ | $0.96$ | $0.00$ | $0.24$ | $0.24$ | $0.91$ | $0.03$ | $0.21$ | $0.22$ | $0.94$ | $0.00$ | $0.21$ | $0.21$ | $0.96$
|  | 0.75 | $0.03$ | $0.24$ | $0.25$ | $0.96$ | $0.00$ | $0.24$ | $0.24$ | $0.91$ | $0.02$ | $0.20$ | $0.20$ | $0.94$ | $0.00$ | $0.20$ | $0.20$ | $0.96$
|  | 1.00 | $0.03$ | $0.24$ | $0.25$ | $0.96$ | $0.00$ | $0.24$ | $0.24$ | $0.91$ | $0.02$ | $0.19$ | $0.19$ | $0.94$ | $0.00$ | $0.19$ | $0.19$ | $0.95$
12 | 12 | 0.50 | $0.09$ | $0.41$ | $0.45$ | $0.96$ | $-0.01$ | $0.41$ | $0.37$ | $0.89$ | $0.08$ | $0.35$ | $0.37$ | $0.94$ | $-0.02$ | $0.35$ | $0.34$ | $0.98$
|  | 0.75 | $0.09$ | $0.41$ | $0.45$ | $0.96$ | $-0.01$ | $0.41$ | $0.37$ | $0.89$ | $0.06$ | $0.33$ | $0.34$ | $0.94$ | $-0.01$ | $0.33$ | $0.32$ | $0.96$
|  | 1.00 | $0.09$ | $0.41$ | $0.45$ | $0.96$ | $-0.01$ | $0.41$ | $0.37$ | $0.89$ | $0.04$ | $0.31$ | $0.31$ | $0.95$ | $-0.01$ | $0.31$ | $0.29$ | $0.97$
| 24 | 0.50 | $0.04$ | $0.28$ | $0.31$ | $0.96$ | $0.00$ | $0.28$ | $0.28$ | $0.93$ | $0.03$ | $0.24$ | $0.24$ | $0.94$ | $-0.01$ | $0.24$ | $0.23$ | $0.96$
|  | 0.75 | $0.04$ | $0.28$ | $0.31$ | $0.96$ | $0.00$ | $0.28$ | $0.28$ | $0.93$ | $0.03$ | $0.23$ | $0.23$ | $0.95$ | $0.00$ | $0.23$ | $0.22$ | $0.96$
|  | 1.00 | $0.04$ | $0.28$ | $0.31$ | $0.96$ | $0.00$ | $0.28$ | $0.28$ | $0.93$ | $0.02$ | $0.21$ | $0.21$ | $0.95$ | $-0.01$ | $0.21$ | $0.21$ | $0.96$
| 32 | 0.50 | $0.02$ | $0.24$ | $0.24$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.92$ | $0.02$ | $0.21$ | $0.21$ | $0.95$ | $-0.01$ | $0.21$ | $0.20$ | $0.95$
|  | 0.75 | $0.02$ | $0.24$ | $0.24$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.92$ | $0.02$ | $0.19$ | $0.20$ | $0.94$ | $-0.01$ | $0.19$ | $0.19$ | $0.95$
|  | 1.00 | $0.02$ | $0.24$ | $0.24$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.92$ | $0.01$ | $0.18$ | $0.18$ | $0.95$ | $0.00$ | $0.18$ | $0.18$ | $0.95$
18 | 12 | 0.50 | $0.06$ | $0.41$ | $0.42$ | $0.97$ | $-0.03$ | $0.41$ | $0.35$ | $0.89$ | $0.07$ | $0.35$ | $0.35$ | $0.95$ | $-0.02$ | $0.35$ | $0.32$ | $0.98$
|  | 0.75 | $0.06$ | $0.41$ | $0.42$ | $0.97$ | $-0.03$ | $0.41$ | $0.35$ | $0.89$ | $0.05$ | $0.32$ | $0.32$ | $0.96$ | $-0.02$ | $0.32$ | $0.30$ | $0.97$
|  | 1.00 | $0.06$ | $0.41$ | $0.42$ | $0.97$ | $-0.03$ | $0.41$ | $0.35$ | $0.89$ | $0.03$ | $0.30$ | $0.30$ | $0.95$ | $-0.01$ | $0.30$ | $0.28$ | $0.98$
| 24 | 0.50 | $0.03$ | $0.28$ | $0.30$ | $0.96$ | $-0.01$ | $0.28$ | $0.28$ | $0.93$ | $0.03$ | $0.24$ | $0.25$ | $0.94$ | $-0.01$ | $0.24$ | $0.24$ | $0.96$
|  | 0.75 | $0.03$ | $0.28$ | $0.30$ | $0.96$ | $-0.01$ | $0.28$ | $0.28$ | $0.93$ | $0.02$ | $0.22$ | $0.23$ | $0.94$ | $-0.01$ | $0.22$ | $0.22$ | $0.96$
|  | 1.00 | $0.03$ | $0.28$ | $0.30$ | $0.96$ | $-0.01$ | $0.28$ | $0.28$ | $0.93$ | $0.02$ | $0.21$ | $0.22$ | $0.94$ | $-0.01$ | $0.21$ | $0.21$ | $0.95$
| 32 | 0.50 | $0.03$ | $0.24$ | $0.25$ | $0.96$ | $0.00$ | $0.24$ | $0.23$ | $0.90$ | $0.03$ | $0.21$ | $0.21$ | $0.94$ | $0.00$ | $0.21$ | $0.20$ | $0.96$
|  | 0.75 | $0.03$ | $0.24$ | $0.25$ | $0.96$ | $0.00$ | $0.24$ | $0.23$ | $0.90$ | $0.02$ | $0.19$ | $0.20$ | $0.94$ | $0.00$ | $0.19$ | $0.19$ | $0.95$
|  | 1.00 | $0.03$ | $0.24$ | $0.25$ | $0.96$ | $0.00$ | $0.24$ | $0.23$ | $0.90$ | $0.02$ | $0.18$ | $0.19$ | $0.94$ | $0.00$ | $0.18$ | $0.19$ | $0.95$

  * •

Note: Bias and ESE are, respectively, the empirical relative bias and standard
error of the IUPM estimator; ASE is the average of the standard error
estimator; CP is the empirical coverage probability of the 95% confidence
interval for the IUPM. There were a total of 39 simulated assays out of
$\mathrm{27{,}000}$ ($0.1\%$) where the MLE and bias-corrected MLE without
UDSA were infinite that were excluded; all other entries are based on
$\mathrm{1000}$ replicates.

Table S2: Simulation results with a single dilution level and a constant rate
of infected cells for all distinct viral lineages. The true overall IUPM in
all settings was $T=0.5$.

| Without UDSA | With UDSA
---|---|---
| MLE | Bias-Corrected MLE | MLE | Bias-Corrected MLE
$\boldsymbol{n^{\prime}}$ | $\boldsymbol{M}$ | $\boldsymbol{q}$ | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP
6 | 12 | 0.50 | $0.03$ | $0.24$ | $0.26$ | $0.95$ | $0.00$ | $0.24$ | $0.23$ | $0.75$ | $0.04$ | $0.23$ | $0.24$ | $0.95$ | $-0.02$ | $0.23$ | $0.22$ | $0.98$
|  | 0.75 | $0.03$ | $0.24$ | $0.26$ | $0.95$ | $0.00$ | $0.24$ | $0.23$ | $0.75$ | $0.03$ | $0.22$ | $0.23$ | $0.96$ | $-0.01$ | $0.22$ | $0.22$ | $0.98$
|  | 1.00 | $0.03$ | $0.24$ | $0.26$ | $0.95$ | $0.00$ | $0.24$ | $0.23$ | $0.75$ | $0.02$ | $0.21$ | $0.22$ | $0.95$ | $0.00$ | $0.21$ | $0.21$ | $0.98$
| 24 | 0.50 | $0.01$ | $0.17$ | $0.17$ | $0.94$ | $0.00$ | $0.17$ | $0.17$ | $0.76$ | $0.01$ | $0.16$ | $0.16$ | $0.95$ | $-0.01$ | $0.16$ | $0.16$ | $0.97$
|  | 0.75 | $0.01$ | $0.17$ | $0.17$ | $0.94$ | $0.00$ | $0.17$ | $0.17$ | $0.76$ | $0.01$ | $0.15$ | $0.16$ | $0.95$ | $-0.01$ | $0.15$ | $0.15$ | $0.97$
|  | 1.00 | $0.01$ | $0.17$ | $0.17$ | $0.94$ | $0.00$ | $0.17$ | $0.17$ | $0.76$ | $0.01$ | $0.15$ | $0.15$ | $0.95$ | $0.00$ | $0.15$ | $0.15$ | $0.96$
| 32 | 0.50 | $0.01$ | $0.14$ | $0.15$ | $0.95$ | $0.00$ | $0.14$ | $0.14$ | $0.69$ | $0.01$ | $0.14$ | $0.14$ | $0.94$ | $-0.01$ | $0.14$ | $0.13$ | $0.95$
|  | 0.75 | $0.01$ | $0.14$ | $0.15$ | $0.95$ | $0.00$ | $0.14$ | $0.14$ | $0.69$ | $0.01$ | $0.13$ | $0.13$ | $0.94$ | $0.00$ | $0.13$ | $0.13$ | $0.95$
|  | 1.00 | $0.01$ | $0.14$ | $0.15$ | $0.95$ | $0.00$ | $0.14$ | $0.14$ | $0.69$ | $0.01$ | $0.13$ | $0.13$ | $0.94$ | $0.00$ | $0.13$ | $0.13$ | $0.95$
12 | 12 | 0.50 | $0.02$ | $0.24$ | $0.25$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.76$ | $0.03$ | $0.23$ | $0.23$ | $0.95$ | $-0.03$ | $0.23$ | $0.21$ | $0.98$
|  | 0.75 | $0.02$ | $0.24$ | $0.25$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.76$ | $0.02$ | $0.22$ | $0.22$ | $0.96$ | $-0.01$ | $0.22$ | $0.21$ | $0.98$
|  | 1.00 | $0.02$ | $0.24$ | $0.25$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.76$ | $0.02$ | $0.21$ | $0.21$ | $0.97$ | $0.00$ | $0.21$ | $0.20$ | $0.97$
| 24 | 0.50 | $0.01$ | $0.17$ | $0.17$ | $0.93$ | $0.00$ | $0.17$ | $0.17$ | $0.73$ | $0.02$ | $0.16$ | $0.16$ | $0.95$ | $-0.01$ | $0.16$ | $0.15$ | $0.96$
|  | 0.75 | $0.01$ | $0.17$ | $0.17$ | $0.93$ | $0.00$ | $0.17$ | $0.17$ | $0.73$ | $0.01$ | $0.15$ | $0.15$ | $0.96$ | $-0.01$ | $0.15$ | $0.15$ | $0.97$
|  | 1.00 | $0.01$ | $0.17$ | $0.17$ | $0.93$ | $0.00$ | $0.17$ | $0.17$ | $0.73$ | $0.01$ | $0.15$ | $0.15$ | $0.96$ | $0.00$ | $0.15$ | $0.14$ | $0.97$
| 32 | 0.50 | $0.01$ | $0.14$ | $0.15$ | $0.95$ | $0.00$ | $0.14$ | $0.14$ | $0.70$ | $0.01$ | $0.14$ | $0.14$ | $0.95$ | $-0.01$ | $0.14$ | $0.13$ | $0.96$
|  | 0.75 | $0.01$ | $0.14$ | $0.15$ | $0.95$ | $0.00$ | $0.14$ | $0.14$ | $0.70$ | $0.01$ | $0.13$ | $0.13$ | $0.95$ | $0.00$ | $0.13$ | $0.13$ | $0.96$
|  | 1.00 | $0.01$ | $0.14$ | $0.15$ | $0.95$ | $0.00$ | $0.14$ | $0.14$ | $0.70$ | $0.01$ | $0.13$ | $0.13$ | $0.96$ | $0.00$ | $0.13$ | $0.12$ | $0.96$
18 | 12 | 0.50 | $0.02$ | $0.24$ | $0.25$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.76$ | $0.03$ | $0.23$ | $0.23$ | $0.95$ | $-0.02$ | $0.23$ | $0.22$ | $0.97$
|  | 0.75 | $0.02$ | $0.24$ | $0.25$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.76$ | $0.02$ | $0.22$ | $0.22$ | $0.96$ | $-0.01$ | $0.22$ | $0.21$ | $0.97$
|  | 1.00 | $0.02$ | $0.24$ | $0.25$ | $0.95$ | $-0.01$ | $0.24$ | $0.23$ | $0.76$ | $0.02$ | $0.21$ | $0.21$ | $0.96$ | $-0.01$ | $0.21$ | $0.20$ | $0.97$
| 24 | 0.50 | $0.01$ | $0.17$ | $0.17$ | $0.94$ | $-0.01$ | $0.17$ | $0.17$ | $0.73$ | $0.01$ | $0.16$ | $0.16$ | $0.95$ | $-0.01$ | $0.16$ | $0.16$ | $0.96$
|  | 0.75 | $0.01$ | $0.17$ | $0.17$ | $0.94$ | $-0.01$ | $0.17$ | $0.17$ | $0.73$ | $0.01$ | $0.15$ | $0.15$ | $0.95$ | $-0.01$ | $0.15$ | $0.15$ | $0.96$
|  | 1.00 | $0.01$ | $0.17$ | $0.17$ | $0.94$ | $-0.01$ | $0.17$ | $0.17$ | $0.73$ | $0.01$ | $0.15$ | $0.15$ | $0.96$ | $0.00$ | $0.15$ | $0.15$ | $0.96$
| 32 | 0.50 | $0.01$ | $0.14$ | $0.15$ | $0.94$ | $0.00$ | $0.14$ | $0.14$ | $0.70$ | $0.01$ | $0.14$ | $0.13$ | $0.94$ | $0.00$ | $0.14$ | $0.13$ | $0.95$
|  | 0.75 | $0.01$ | $0.14$ | $0.15$ | $0.94$ | $0.00$ | $0.14$ | $0.14$ | $0.70$ | $0.01$ | $0.13$ | $0.13$ | $0.94$ | $0.00$ | $0.13$ | $0.13$ | $0.96$
|  | 1.00 | $0.01$ | $0.14$ | $0.15$ | $0.94$ | $0.00$ | $0.14$ | $0.14$ | $0.70$ | $0.01$ | $0.13$ | $0.13$ | $0.96$ | $0.00$ | $0.13$ | $0.13$ | $0.96$

  * •

Note: Bias and ESE are, respectively, the empirical relative bias and standard
error of the IUPM estimator; ASE is the average of the standard error
estimator; CP is the empirical coverage probability of the 95% confidence
interval for the IUPM. All entries are based on $\mathrm{1000}$ replicates.

Table S3: Simulation results with multiple dilution levels and a non-constant
rate of infected cells for all distinct viral lineages. In all settings, the
true IUPM was $T=1$ and the proportions of p24-positive wells that were deep-
sequenced at the three dilution levels were $\boldsymbol{q}=(0,0.5,1)$.

| Without UDSA | With UDSA
---|---|---
| MLE | Bias-Corrected MLE | MLE | Bias-Corrected MLE
$\boldsymbol{n^{\prime}}$ | $\boldsymbol{M}$ | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP | Bias | ASE | ESE | CP
6 | 6, 12, 18 | $0.04$ | $0.22$ | $0.24$ | $0.95$ | $0.00$ | $0.22$ | $0.22$ | $0.96$ | $0.02$ | $0.16$ | $0.16$ | $0.94$ | $-0.01$ | $0.16$ | $0.16$ | $0.95$
| 9, 18, 27 | $0.03$ | $0.18$ | $0.19$ | $0.95$ | $0.01$ | $0.18$ | $0.18$ | $0.96$ | $0.02$ | $0.13$ | $0.13$ | $0.96$ | $0.00$ | $0.13$ | $0.13$ | $0.96$
| 12, 24, 36 | $0.02$ | $0.16$ | $0.16$ | $0.95$ | $0.00$ | $0.16$ | $0.16$ | $0.96$ | $0.01$ | $0.11$ | $0.11$ | $0.96$ | $0.00$ | $0.11$ | $0.11$ | $0.96$
12 | 6, 12, 18 | $0.03$ | $0.22$ | $0.24$ | $0.96$ | $-0.01$ | $0.22$ | $0.22$ | $0.96$ | $0.02$ | $0.15$ | $0.16$ | $0.95$ | $0.00$ | $0.15$ | $0.15$ | $0.95$
| 9, 18, 27 | $0.03$ | $0.18$ | $0.18$ | $0.96$ | $0.00$ | $0.18$ | $0.18$ | $0.96$ | $0.02$ | $0.13$ | $0.13$ | $0.94$ | $0.00$ | $0.13$ | $0.12$ | $0.95$
| 12, 24, 36 | $0.02$ | $0.16$ | $0.16$ | $0.95$ | $0.00$ | $0.16$ | $0.16$ | $0.96$ | $0.01$ | $0.11$ | $0.11$ | $0.95$ | $0.00$ | $0.11$ | $0.11$ | $0.95$
18 | 9, 18, 27 | $0.02$ | $0.18$ | $0.19$ | $0.95$ | $0.00$ | $0.18$ | $0.18$ | $0.95$ | $0.01$ | $0.12$ | $0.12$ | $0.95$ | $0.00$ | $0.12$ | $0.12$ | $0.95$
| 6, 12, 18 | $0.04$ | $0.22$ | $0.24$ | $0.96$ | $0.00$ | $0.22$ | $0.22$ | $0.96$ | $0.02$ | $0.15$ | $0.16$ | $0.94$ | $0.00$ | $0.15$ | $0.15$ | $0.95$
| 12, 24, 36 | $0.02$ | $0.16$ | $0.16$ | $0.96$ | $0.00$ | $0.16$ | $0.15$ | $0.96$ | $0.01$ | $0.11$ | $0.11$ | $0.94$ | $0.00$ | $0.11$ | $0.11$ | $0.95$

  * •

Note: Bias and ESE are, respectively, the empirical relative bias and standard
error of the IUPM estimator; ASE is the average of the standard error
estimator; CP is the empirical coverage probability of the 95% confidence
interval for the IUPM. All entries are based on $\mathrm{1000}$ replicates.

Table S4: Summary of assay data collected from 17 individuals receiving care
at the University of North Carolina HIV Cure Center

| QVOA | UDSA
---|---|---
ID | $\boldsymbol{u}$ | $\boldsymbol{M}$ | $\boldsymbol{M_{P}}$ | $\boldsymbol{m}$ | $\boldsymbol{n}$
C1 | $(2.5,0.5,0.1)$ | $(36,6,6)$ | $(4,0,0)$ | $(4,0,0)$ | 12
C2 | $(2.5,0.5,0.1)$ | $(36,6,6)$ | $(5,1,0)$ | $(5,0,0)$ | 7
C3 | $(2.5,0.5,0.1)$ | $(18,6,6)$ | $(5,0,0)$ | $(5,0,0)$ | 4
C4 | $(2.5,0.5,0.1)$ | $(18,6,6)$ | $(5,1,0)$ | $(3,0,0)$ | 3
C5 | $(2.5,0.5,0.1,0.025)$ | $(14,6,6,6)$ | $(4,0,1,0)$ | $(3,0,0,0)$ | 3
C6 | $(2.5,0.5,0.1)$ | $(18,6,6)$ | $(7,0,0)$ | $(6,0,0)$ | 4
C7 | $(2.5,0.5,0.1)$ | $(36,6,6)$ | $(15,1,0)$ | $(15,0,0)$ | 39
C8 | $(2.5,0.5,0.1)$ | $(36,6,6)$ | $(22,2,1)$ | $(22,0,0)$ | 26
C9 | $(2.5,0.5,0.1,0.025)$ | $(12,6,6,6)$ | $(9,0,0,0)$ | $(6,0,0,0)$ | 8
C10 | $(2.5,0.5,0.1)$ | $(18,6,6)$ | $(12,3,1)$ | $(6,0,0)$ | 19
C11 | $(2.5,0.5,0.1,0.025)$ | $(12,6,6,6)$ | $(9,1,1,1)$ | $(6,0,0,0)$ | 9
C12 | $(2.5,0.5,0.1)$ | $(36,6,6)$ | $(32,3,1)$ | $(32,0,0)$ | 65
C13 | $(2.5,0.5,0.1,0.025)$ | $(18,6,6,6)$ | $(16,4,3,0)$ | $(0,4,0,0)$ | 7
C14 | $(2.5,0.5,0.1)$ | $(18,6,6)$ | $(18,3,1)$ | $(0,3,0)$ | 3
C15 | $(2.5,0.5,0.1)$ | $(18,6,6)$ | $(18,5,0)$ | $(0,5,0)$ | 6
C16 | $(2.5,0.5,0.1,0.025)$ | $(12,6,6,6)$ | $(12,4,2,0)$ | $(0,4,0,0)$ | 8
C17 | $(2.5,0.5,0.1,0.025)$ | $(18,6,6,6)$ | $(18,4,3,1)$ | $(0,4,0,0)$ | 6

  * •

Note: $\boldsymbol{u}$ is the vector of dilution levels in millions of cells
per well; $\boldsymbol{M}$ is the vector of replicate well numbers to undergo
the QVOA at each dilution level; $\boldsymbol{M_{P}}$ is the vector of
positive wells found by the QVOA at each dilution level; $\boldsymbol{m}$ is
the vector of positive well numbers to undergo the UDSA at each dilution
level; $\boldsymbol{n}$ is the number of distinct viral lineages (DVL)
detected by the UDSA.

Table S5: Estimated infectious units per million (IUPM) of HIV (with 95% confidence intervals) for 17 individuals receiving care at the University of North Carolina HIV Cure Center | Method
---|---
Subject | Bias | Without UDSA | With UDSA | With UDSA
ID | Correction | (Multiple Dilutions) | (Single Dilution) | (Multiple Dilutions)
C1 | MLE |  $0.05$ $(0.02,0.12)$ |  $0.14$ $(0.08,0.24)$ |  $0.13$ $(0.07,0.23)$
C1 | BC-MLE |  $0.04$ $(0.02,0.12)$ |  $0.13$ $(0.08,0.24)$ |  $0.13$ $(0.07,0.23)$
C2 | MLE |  $0.07$ $(0.03,0.15)$ |  $0.17$ $(0.10,0.29)$ |  $0.18$ $(0.11,0.29)$
C2 | BC-MLE |  $0.07$ $(0.03,0.15)$ |  $0.17$ $(0.10,0.29)$ |  $0.18$ $(0.11,0.29)$
C3 | MLE |  $0.12$ $(0.05,0.29)$ |  $0.12$ $(0.05,0.28)$ |  $0.11$ $(0.04,0.26)$
C3 | BC-MLE |  $0.12$ $(0.05,0.29)$ |  $0.11$ $(0.05,0.28)$ |  $0.10$ $(0.04,0.26)$
C4 | MLE |  $0.14$ $(0.06,0.32)$ |  $0.12$ $(0.05,0.29)$ |  $0.14$ $(0.06,0.30)$
C4 | BC-MLE |  $0.14$ $(0.06,0.32)$ |  $0.11$ $(0.04,0.29)$ |  $0.13$ $(0.05,0.30)$
C5 | MLE |  $0.15$ $(0.06,0.36)$ |  $0.20$ $(0.09,0.44)$ |  $0.20$ $(0.10,0.43)$
C5 | BC-MLE |  $0.14$ $(0.06,0.36)$ |  $0.18$ $(0.08,0.44)$ |  $0.19$ $(0.08,0.43)$
C6 | MLE |  $0.18$ $(0.08,0.38)$ |  $0.34$ $(0.20,0.60)$ |  $0.31$ $(0.18,0.54)$
C6 | BC-MLE |  $0.17$ $(0.08,0.37)$ |  $0.33$ $(0.19,0.59)$ |  $0.30$ $(0.17,0.53)$
C7 | MLE |  $0.22$ $(0.13,0.36)$ |  $0.64$ $(0.49,0.83)$ |  $0.63$ $(0.48,0.81)$
C7 | BC-MLE |  $0.22$ $(0.13,0.36)$ |  $0.63$ $(0.48,0.82)$ |  $0.62$ $(0.48,0.81)$
C8 | MLE |  $0.41$ $(0.27,0.62)$ |  $1.06$ $(0.86,1.31)$ |  $1.06$ $(0.86,1.30)$
C8 | BC-MLE |  $0.40$ $(0.26,0.61)$ |  $1.04$ $(0.84,1.29)$ |  $1.04$ $(0.84,1.29)$
C9 | MLE |  $0.44$ $(0.22,0.87)$ |  $0.60$ $(0.35,1.03)$ |  $0.51$ $(0.30,0.88)$
C9 | BC-MLE |  $0.42$ $(0.20,0.86)$ |  $0.56$ $(0.31,1.00)$ |  $0.48$ $(0.28,0.85)$
C10 | MLE |  $0.54$ $(0.32,0.91)$ |  $0.73$ $(0.48,1.12)$ |  $0.79$ $(0.53,1.17)$
C10 | BC-MLE |  $0.52$ $(0.30,0.89)$ |  $0.69$ $(0.44,1.09)$ |  $0.75$ $(0.50,1.14)$
C11 | MLE |  $0.62$ $(0.34,1.14)$ |  $0.89$ $(0.55,1.43)$ |  $0.88$ $(0.57,1.38)$
C11 | BC-MLE |  $0.58$ $(0.30,1.11)$ |  $0.82$ $(0.49,1.37)$ |  $0.83$ $(0.51,1.33)$
C12 | MLE |  $0.94$ $(0.63,1.39)$ |  $1.42$ $(1.19,1.69)$ |  $1.42$ $(1.19,1.69)$
C12 | BC-MLE |  $0.90$ $(0.59,1.36)$ |  $1.40$ $(1.17,1.67)$ |  $1.40$ $(1.17,1.67)$
C13 | MLE |  $1.24$ $(0.74,2.08)$ |  $2.55$ $(1.22,5.36)$ |  $1.41$ $(0.79,2.49)$
C13 | BC-MLE |  $1.12$ $(0.63,1.99)$ |  $2.32$ $(1.03,5.25)$ |  $1.19$ $(0.61,2.34)$
C14 | MLE |  $1.82$ $(0.90,3.66)$ |  $2.12$ $(0.86,5.18)$ |  $2.20$ $(1.12,4.34)$
C14 | BC-MLE |  $1.54$ $(0.67,3.52)$ |  $1.88$ $(0.69,5.15)$ |  $1.73$ $(0.73,4.11)$
C15 | MLE |  $2.49$ $(1.10,5.62)$ |  $3.66$ $(1.88,7.10)$ |  $2.95$ $(1.57,5.57)$
C15 | BC-MLE |  $2.08$ $(0.79,5.52)$ |  $3.27$ $(1.56,6.87)$ |  $2.62$ $(1.28,5.36)$
C16 | MLE |  $2.52$ $(1.16,5.46)$ |  $2.92$ $(1.46,5.84)$ |  $2.97$ $(1.62,5.44)$
C16 | BC-MLE |  $2.18$ $(0.89,5.33)$ |  $2.65$ $(1.24,5.69)$ |  $2.70$ $(1.39,5.26)$
C17 | MLE |  $3.49$ $(1.69,7.20)$ |  $2.19$ $(0.98,4.88)$ |  $3.03$ $(1.66,5.55)$
C17 | BC-MLE |  $3.08$ $(1.36,7.00)$ |  $1.99$ $(0.82,4.80)$ |  $2.74$ $(1.40,5.35)$

##  Web Appendix E Derivations and Simulations for the Negative Binomial MLE

###  Web Appendix E.1 Negative binomial likelihood

Suppose $X_{ij}$ follows a negative binomial distribution with mean parameter
$\lambda_{i}$ and dispersion parameter $\gamma$ (shared across all DVLs),
i.e., $X_{ij}$ has probability mass function (PMF)

$\displaystyle{\Pr}_{\lambda_{i},\gamma}(X_{ij}=x_{ij})=\frac{\Gamma(x_{ij}+1/\gamma)}{\Gamma(1/\gamma)\Gamma(x_{ij}+1)}\left(\frac{1}{\gamma\lambda_{i}+1}\right)^{1/\gamma}\left(\frac{\gamma\lambda_{i}}{\gamma\lambda_{i}+1}\right)^{x_{ij}},$

for $\gamma\in(0,\infty)$, $x_{ij}\in\\{0,1,2,...\\}$, and $\Gamma(\cdot)$
denoting the gamma function (agresti_categorical_2013). To allow
$\gamma\in[0,\infty)$, define ${\Pr}_{\lambda_{i},\gamma=0}(X_{ij}=x_{ij})$ as

$\displaystyle\lim_{\gamma\to
0}{\Pr}_{\lambda_{i},\gamma}(X_{ij}=x_{ij})=\frac{\exp(-\lambda_{i})\lambda_{i}^{x_{ij}}}{x_{ij}!}.$

That is, if $\gamma=0$ then $X_{ij}$ follows a Poisson distribution with mean
$\lambda_{i}$, as in Section 2 of the main text. When $\gamma>0$,
${\Pr}_{\lambda_{i},\gamma}(X_{ij}=0)=(\gamma\lambda_{i}+1)^{-1/\gamma}$, so
$Z_{ij}$ and $W_{j}$ follow Bernoulli distributions with
$\Pr_{\lambda_{i},\gamma}(Z_{ij}=1)=1-(\gamma\lambda_{i}+1)^{-1/\gamma}$ and
$\Pr_{\boldsymbol{\lambda},\gamma}(W_{j}=1)=1-\prod_{i=1}^{n}(\gamma\lambda_{i}+1)^{-1/\gamma}$.
Following Sections 2.2 and 2.6 of the main text, when $\gamma>0$, the negative
binomial observed-data likelihood for data from a single dilution level is
proportional to $L_{NB}(\boldsymbol{\lambda},\gamma|M_{N},\boldsymbol{Y})=$

$\displaystyle\prod_{i=1}^{n}\left[\left\\{1-(\gamma\lambda_{i}+1)^{-1/\gamma}\right\\}^{Y_{i}}(\gamma\lambda_{i}+1)^{-(M_{N}+m-Y_{i})/\gamma}\right]\left\\{1-\prod_{i=1}^{n}(\gamma\lambda_{i}+1)^{-1/\gamma}\right\\}^{M-M_{N}-m},$

and for data from multiple dilution levels the likelihood is proportional to

$\displaystyle\widetilde{L}_{NB}(\boldsymbol{\tau},\gamma|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})=\prod_{d=1}^{D}L_{NB}\left(u^{(d)}\boldsymbol{\tau},\gamma|M_{N}^{(d)},\boldsymbol{Y}^{(d)}\right).$
(S.16)

To allow $\gamma\in[0,\infty)$, define
$L_{NB}(\boldsymbol{\lambda},\gamma=0|M_{N},\boldsymbol{Y})\equiv\lim_{\gamma\to
0}L_{NB}(\boldsymbol{\lambda},\gamma|M_{N},\boldsymbol{Y})=L(\boldsymbol{\lambda}|M_{N},\boldsymbol{Y})$
for the Poisson distribution in (1) of the main text.

Note that the negative binomial likelihood for a single dilution level $u$
depends on the parameters $\lambda_{i}=u\tau_{i}$ and $\gamma$ only through
$\Pr_{\lambda_{i},\gamma}(Z_{ij}=1)$, and that different values of
$\lambda_{i}$ and $\gamma$ can lead to the same
$\Pr_{\lambda_{i},\gamma}(Z_{ij}=1)$. For example,
$\Pr_{\lambda_{i}=1,\gamma=2}(Z_{ij}=1)=\Pr_{\lambda_{i}=\sqrt{3}-1,\gamma=1}(Z_{ij}=1)=1-1/\sqrt{3}$.
Thus, the negative binomial parameters are not identifiable given assay data
from a single dilution level. On the other hand, it is straightforward to show
that the negative binomial parameters are identifiable given assay data from
multiple dilution levels. In this case, the negative binomial MLE
($\widehat{\boldsymbol{\tau}}_{NB},\widehat{\gamma}_{NB})$ is obtained by
maximizing the log of (S.16) over $\boldsymbol{\tau}\in(0,\infty)^{n}$ and
$\gamma\in[0,\infty)$ using the L-BFGS-B method (ByrdEtAl1995), a variation of
BFGS that allows for the inequality constraint $\gamma\geq 0$. ByrdEtAl1995
suggest that the L-BFGS-B method will typically converge to a global maximum
given a concave objective function.

###  Web Appendix E.2 Likelihood ratio test for overdispersion

The presence of overdispersion can be assessed by testing the null hypothesis
$H_{0}:\gamma=0$ versus the alternative $H_{1}:\gamma>0$. These hypotheses can
be tested with a likelihood ratio test (LRT) using the test statistic
$S_{LR}=2\log\\{\widetilde{L}_{NB}(\hat{\boldsymbol{\tau}}_{NB},\widehat{\gamma}_{NB}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})/\widetilde{L}(\hat{\boldsymbol{\tau}}|\boldsymbol{M_{N}},\boldsymbol{Y},\boldsymbol{u})\\}$,
comparing the maximum values of the Poisson and negative binomial likelihood
functions. Since this is a test of one parameter ($\gamma$) with the null
value ($0$) on the boundary of the parameter space, under $H_{0}$ the LRT
statistic $S_{LR}$ has an asymptotic distribution that is a 50/50 mixture of
$\chi^{2}_{0}$ and $\chi^{2}_{1}$ distributions (self_asymptotic_1987).

###  Web Appendix E.3 Overdispersion simulations

The performance of the negative binomial MLE and the LRT test for
overdispersion was investigated in a simulation study. Assay data were
simulated following Section 3.2 for $D=3$ dilution levels. The following
parameters were held fixed: (i) the true IUPM $T=1$, (ii) the three dilution
levels $\boldsymbol{u}=(u_{1},u_{2},u_{3})=(0.5,1,2)$ million cells per well,
(iii) the proportions of positive wells to be deep sequenced at the three
dilution levels $\boldsymbol{q}=(q_{1},q_{2},q_{3})=(0,0.5,1)$, and (iv)
$n^{\prime}=12$ underlying DVLs. The simulation settings varied by the number
of replicate wells per dilution level,
$\boldsymbol{M}=(M_{1},M_{2},M_{3})=(6,12,18)$ or $(30,60,90)$, and the
dispersion parameter $\gamma=0$ (no overdispersion), $\gamma=0.2$ (light
overdispersion), $\gamma=1$ (moderate overdispersion), and $\gamma=4$ (heavy
overdispersion). For each scenario, $500$ assays were simulated.

Figure 1 shows the empirical distributions of three IUPM estimators under
varying degrees of overdispersion: (i) the negative binomial MLE, (ii) the
Poisson MLE, and (iii) the Poisson BC-MLE. As expected, when the counts of
infected cells $X_{ij}$ follow a negative binomial distribution with nonzero
$\gamma$, the Poisson MLE and Poisson BC-MLE are biased, regardless of the
numbers of replicate wells $\boldsymbol{M}$. For light to moderate
overdispersion, this bias is relatively small, but it becomes severe under
heavy overdispersion. The negative binomial MLE has notably smaller bias under
overdispersion. However, in all scenarios, the variance of the negative
binomial MLE is much larger than the variances of the Poisson MLE and BC-MLE.
In terms of mean-squared error (MSE), a combined measure of bias and variance,
the Poisson MLE and BC-MLE outperform the negative binomial MLE.

Figure 1: Empirical distributions of the estimated IUPM with $D=3$ dilution
levels (0.5, 1, and 2 million cells per well), proportions of positive wells
that were deep sequenced at the three dilution levels of
$\boldsymbol{q}=(0,0.5,1)$, $n^{\prime}=12$ DVLs, IUPM of $T=1$, and
$\mathrm{500}$ simulations per setting. *Number of values beyond the range of
the plot.

Table S6 shows the empirical power of the LRT of overdispersion at the 0.05
significance level against various alternatives. Under the null hypothesis,
the LRT had approximately correct type 1 error rate ($\approx 0.05$) for both
small and large numbers of replicate wells $\boldsymbol{M}$. Under the
alternative hypotheses, the power of the LRT increased as $\gamma$ increased
(i.e., it became farther from the null value of $0$) and as the numbers of
replicate wells $\boldsymbol{M}$ increased. However, in all cases, the power
of the LRT against overdispersion was very low.

Table S6: Empirical power of the likelihood ratio test for overdispersion

M | $\boldsymbol{\gamma}$ | Power of LRT
---|---|---
(6, 12, 18) | 0.0 | 0.046
(6, 12, 18) | 0.2 | 0.064
(6, 12, 18) | 1.0 | 0.098
(6, 12, 18) | 4.0 | 0.104
(30, 60, 90) | 0.0 | 0.042
(30, 60, 90) | 0.2 | 0.042
(30, 60, 90) | 1.0 | 0.082
(30, 60, 90) | 4.0 | 0.216

  * •

Note: $\gamma$ is the negative binomial dispersion parameter; $\gamma=0$
corresponds to the null hypothesis of no overdispersion; Power is the
proportion of simulated replicates where the null hypothesis of no
overdispersion was rejected at the 0.05 significance level.

In summary, these simulation results demonstrate for the scenarios considered
that it is generally difficult to detect overdispersion in the latent counts
of infected cells given the observed well-level infection indicator variables.
Moreover, despite correctly modeling overdispersion, the negative binomial MLE
for the IUPM has much greater MSE than the Poisson MLE, regardless of the true
dispersion parameter, due to its increased variance.

###  Web Appendix E.4 Assessing overdispersion in the HIV application

The HIV data from Section 4 of the main text were assessed for overdispersion
using the likelihood ratio test described in Web Appendix E.2. Table S7 shows
the overdispersion LRT statistics and corresponding p-values, as well as the
Poisson and negative binomial MLEs for each of the 17 participants.

Table S7: Overdispersion LRT results for HIV Application | Estimated IUPM
---|---
Subject ID | LRT Statistic | LRT P-Value | Poisson MLE | Negative Binomial MLE
C1 | 0.000 | 1.000 | 0.130 | 0.130
C2 | 0.174 | 0.338 | 0.178 | 0.315
C3 | 0.000 | 1.000 | 0.107 | 0.107
C4 | 0.412 | 0.261 | 0.135 | 0.350
C5 | 0.492 | 0.241 | 0.203 | 0.862
C6 | 0.000 | 1.000 | 0.314 | 0.314
C7 | 0.000 | 1.000 | 0.627 | 0.627
C8 | 0.000 | 1.000 | 1.057 | 1.057
C9 | 0.000 | 1.000 | 0.515 | 0.515
C10 | 1.439 | 0.115 | 0.791 | 1.971
C11 | 0.043 | 0.418 | 0.885 | 1.052
C12 | 0.010 | 0.460 | 1.418 | 1.510
C13 | 7.469 | 0.003 | 1.406 | 6.176
C14 | 0.000 | 1.000 | 2.203 | 2.203
C15 | 0.000 | 1.000 | 2.953 | 2.953
C16 | 0.000 | 1.000 | 2.968 | 2.968
C17 | 0.483 | 0.244 | 3.034 | 4.679

Aside from participant C13, the LRT tests provide no evidence of
overdispersion and IUPM estimates are similar for the negative binomial and
Poisson models. For participant C13 the LRT p-value is 0.003, suggesting
possible overdispersion. However, this result should be interpreted with
caution as the reported p-values are not adjusted for multiple testing. For
example, if the LRT p-values for all subjects were compared to a 0.05
significance threshold, then the family-wise error rate would be
$1-(1-0.05)^{17}\approx 0.58$. In addition, the results for subject C13 are
sensitive to slight changes in the observed data, as demonstrated in Table S8.

Table S8: Sensitivity Analysis for Subject C13

$\boldsymbol{M_{P}}$ | LRT Statistic | LRT P-Value | Negative Binomial MLE
---|---|---|---
(16, 4, 3, 0) | 7.469 | 0.003 | 6.179
(17, 4, 3, 0) | 4.450 | 0.017 | 5.190
(18, 4, 3, 0) | 0.057 | 0.406 | 3.361
(16, 4, 2, 0) | 5.058 | 0.012 | 4.204
(16, 4, 1, 0) | 3.145 | 0.038 | 2.937

  * •

Note: The true counts of QVOA-positive wells for participant C13 is
$\boldsymbol{M_{P}}$ = (16, 4, 3, 0), which corresponds to the first row.
Subsequent rows show the effect of changing one component of
$\boldsymbol{M_{P}}$ (underlined) by a small amount.

##  Web Appendix F Extension for Imperfect Assay Sensitivity and Specificity

###  Web Appendix F.1 Single dilution likelihood for imperfect assay
sensitivity and specificity

For simplicity, focus on assay data from a single dilution level of one
million cells per well. Let $W_{j}^{*}$ and $Z_{ij}^{*}$ denote error-prone
versions of the true infection indicators $W_{j}$ and $Z_{ij}$, respectively,
obtained using imperfect assays. The true infection indicators $W_{j}$ and
$\boldsymbol{Z}_{j}$ are missing for all wells (i.e., they are latent
variables). Instead, suppose that the imperfect QVOA results
$\boldsymbol{W}^{*}=(W_{1}^{*},\dots,W_{M}^{*})^{\textrm{T}}$ are fully
observed, while the imperfect UDSA results
$\boldsymbol{Z}^{*}=(\boldsymbol{Z}_{1}^{*},\dots,\boldsymbol{Z}_{M}^{*})^{\textrm{T}}$
will have missing data for any unsequenced wells. Let
$M_{P}^{*}=\sum_{j=1}^{M}W_{j}^{*}$ and $M_{N}^{*}=M-M_{P}^{*}$ denote the
numbers of wells found to be HIV-positive and HIV-negative, respectively,
using the imperfect QVOA, and let $m^{*}=\lfloor qM_{P}^{*}\rceil$ denote the
number of positive wells to be deep-sequenced via the imperfect UDSA.
Importantly, $W_{j}^{*}=0$ no longer implies that
$\boldsymbol{Z}_{j}^{*}=\boldsymbol{0}$, so the imperfect UDSA results are
missing for all ($M-m^{*}$) unsequenced wells.

To incorporate latent $W_{j}$ and $\boldsymbol{Z}_{j}$ in addition to
partially missing $\boldsymbol{Z}_{j}^{*}$, begin with a “complete”
observation,
$(R_{j},W_{j}^{*},\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})$, assumed
to be generated from the joint distribution

$\displaystyle\Pr(R_{j},W_{j}^{*},\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})$
$\displaystyle=\Pr(R_{j}|W_{j}^{*},\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})\Pr(W_{j}^{*}|\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j}^{*}|W_{j},\boldsymbol{Z}_{j})\Pr(W_{j}|\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j})$
$\displaystyle=\Pr(R_{j}|W_{j}^{*})\Pr(W_{j}^{*}|\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j}^{*}|W_{j},\boldsymbol{Z}_{j})\Pr(W_{j}|\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j}),$
(S.17)

where $\Pr(R_{j}|W_{j}^{*})=m^{*}/M_{P}^{*}$ by design;
$\Pr(W_{j}^{*}|\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})$ and
$\Pr(\boldsymbol{Z}_{j}^{*}|W_{j},\boldsymbol{Z}_{j})$ are the error
mechanisms of the imperfect QVOA and UDSA, respectively; and
$\Pr(W_{j}|\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j})=\Pr(W_{j},\boldsymbol{Z}_{j})$
is the joint distribution of the true infection statuses from the perfect
assays (discussed in Sections 2.1–2.6). Since the true QVOA result is fully
determined by the true UDSA results,
$\Pr(W_{j}|\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j})=\Pr(\boldsymbol{Z}_{j})$
still.

Further assumptions are made about the error mechanisms. First, the error-
prone assay results are assumed to depend only on the true assay results for
the same test, such that the error mechanisms simplify to
$\Pr(W_{j}^{*}|\boldsymbol{Z}_{j}^{*},W_{j},\boldsymbol{Z}_{j})=\Pr(W_{j}^{*}|W_{j})$
and
$\Pr(\boldsymbol{Z}_{j}^{*}|W_{j},\boldsymbol{Z}_{j})=\Pr(\boldsymbol{Z}_{j}^{*}|\boldsymbol{Z}_{j})$.
Note that $\Pr(W_{j}^{*}=1|W_{j}=1)$ and $\Pr(Z_{ij}^{*}=1|Z_{ij}=1)$ are the
sensitivities and $\Pr(W_{j}^{*}=0|W_{j}=0)$ and $\Pr(Z_{ij}^{*}=0|Z_{ij}=0)$
are specificites of the QVOA and UDSA assays, respectively. Second, the error-
prone UDSA results for the $i$th DVL are assumed to depend only on the true
UDSA results for the same DVL, such that
$\Pr(Z_{ij}^{*}|\boldsymbol{Z}_{j})=\Pr(Z_{ij}^{*}|Z_{ij})$ and thus
$\Pr(\boldsymbol{Z}_{j}^{*}|\boldsymbol{Z}_{j})=\prod_{i=1}^{n}\Pr(Z_{ij}^{*}|Z_{ij}).$
Third, the sensitivity and specificity of both assays are assumed to be known
(e.g., from the literature or a validation study). If sensitivity and
specificity values are not known, then analyses may be repeated over a range
of plausible values to assess robustness of the IUPM estimates against
imperfect assays.

With these assumptions, the joint distribution simplifies from (S.17) to

$\displaystyle\Pr(R_{j}|W_{j}^{*})\Pr(W_{j}^{*}|W_{j})\Pr(\boldsymbol{Z}_{j}^{*}|\boldsymbol{Z}_{j})\Pr(\boldsymbol{Z}_{j}),$
(S.18)

where $\Pr(\boldsymbol{Z}_{j})=\prod_{i=1}^{n}\Pr(Z_{ij})$ following from the
assumption of independence between the counts of infected cells across DVLs,
as in the main text.

To incorporate all available information on all wells, the observed-data
likelihood function for imperfect assays can be constructed from (S.18), and
is proportional to

$\displaystyle
L(\boldsymbol{\lambda}|\boldsymbol{W}^{*},\boldsymbol{Z}^{*},\boldsymbol{R})$
$\displaystyle=\prod_{j=1}^{M}{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*},\boldsymbol{Z}_{j}^{*})^{R_{j}}{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*})^{1-R_{j}},$
(S.19)

where $\Pr_{\boldsymbol{\lambda}}(W_{j}^{*},\boldsymbol{Z}_{j}^{*})$ is the
PMF of the imperfect assay results $(W_{j}^{*},\boldsymbol{Z}_{j}^{*})$ and
$\Pr_{\boldsymbol{\lambda}}(W_{j}^{*})$ is the marginal PMF of the imperfect
QVOA result $W_{j}^{*}$. Notice that (S.19) follows the same general form as
$L(\boldsymbol{\lambda}|\boldsymbol{W},\boldsymbol{Z},\boldsymbol{R})$ under
perfect sensitivity and specificity from the main text, except that now
$\boldsymbol{W}^{*}$ are fully observed (instead of $\boldsymbol{W}$), while
$\boldsymbol{Z}^{*}$ has some missing data (instead of $\boldsymbol{Z}$). The
imperfect UDSA results $\boldsymbol{Z}^{*}$ are still MAR for the $M-m^{*}$
unsequenced wells, so $\Pr(R_{j}|W_{j}^{*})$ can be omitted from the
likelihood for $\boldsymbol{\lambda}$ (Little&Rubin2002).

The PMFs for the imperfect assay results depend on the true assay results and
can be calculated by marginalizing out the unobserved variables to obtain

$\displaystyle{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*},\boldsymbol{Z}_{j}^{*})$
$\displaystyle=\sum_{z_{1}=0}^{1}\cdots\sum_{z_{n}=0}^{1}\prod_{i=1}^{n}\Pr(W_{j}^{*}|W_{j}=w)\Pr(Z_{ij}^{*}|Z_{ij}=z_{i}){\Pr}_{\lambda_{i}}(Z_{ij}=z_{i}),$
(S.20)

for $w=\textrm{I}(\sum_{i=1}^{n}z_{i}\geq 1)$, and

$\displaystyle{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*})$
$\displaystyle=\sum_{w=0}^{1}\Pr(W_{j}^{*}|W_{j}=w){\Pr}_{\boldsymbol{\lambda}}(W_{j}=w),$
(S.21)

where, as in the main text,
${\Pr}_{\lambda_{i}}(Z_{ij})=\left\\{1-\exp(-\lambda_{i})\right\\}^{Z_{ij}}\exp(-\lambda_{i})^{(1-Z_{ij})}$
and
${\Pr}_{\boldsymbol{\lambda}}(W_{j})=\left\\{1-\exp\left(-\sum_{i=1}^{n}\lambda_{i}\right)\right\\}^{W_{j}}\exp\left(-\sum_{i=1}^{n}\lambda_{i}\right)^{1-W_{j}}$.
Note that when there is no assay error, (S.20) and (S.21) simplify to
${\Pr}_{\boldsymbol{\lambda}}(W_{j},\boldsymbol{Z}_{j})$ and
${\Pr}_{\boldsymbol{\lambda}}(W_{j})$, respectively, from Section 2.2.

When the assays have imperfect sensitivity and specificity, the MLE
$\widehat{\lambda}_{i}$ for detected DVL $i$ is not guaranteed to be positive.
In this case, the L-BFGS-B method introduced in Web Appendix E is used to
impose the inequality constraint $\lambda_{i}\geq 0$. As in Web Appendix A,
supplying the gradient of the log-likelihood based on (S.19) improves the
performance of the numerical optimization procedure. The elements of this
gradient are given by

$\displaystyle\frac{\partial}{\partial\lambda_{i}}l(\boldsymbol{\lambda}|\boldsymbol{W}^{*},\boldsymbol{Z}^{*},\boldsymbol{R})$
$\displaystyle=\sum_{j=1}^{M}R_{j}\frac{\frac{\partial}{\partial\lambda_{i}}{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*},\boldsymbol{Z}_{j}^{*})}{{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*},\boldsymbol{Z}_{j}^{*})}+\sum_{j=1}^{M}(1-R_{j})\frac{{\frac{\partial}{\partial\lambda_{i}}\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*})}{{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*})},$

where
$\frac{\partial}{\partial\lambda_{i}}{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*},\boldsymbol{Z}_{j}^{*})=$

$\displaystyle\sum_{z_{1}=0}^{1}\cdots\sum_{z_{n}=0}^{1}\Pr(W_{j}^{*}|W_{j}=w)\Pr(\boldsymbol{Z}_{j}^{*}|\boldsymbol{Z}_{j}=\boldsymbol{z})\left\\{\prod_{i^{\prime}\neq
i}{\Pr}_{\lambda_{i^{\prime}}}(Z_{i^{\prime}j}=z_{i^{\prime}})\right\\}\left\\{(-1)^{1-z_{i}}\exp(-\lambda_{i})\right\\}$

and
$\frac{\partial}{\partial\lambda_{i}}{\Pr}_{\boldsymbol{\lambda}}(W_{j}^{*})=$

$\displaystyle\sum_{w=0}^{1}\Pr(W_{j}^{*}|W_{j}=w)(-1)^{1-w}\exp\left(-\sum_{i^{\prime}=1}^{n}\lambda_{i^{\prime}}\right).$

###  Web Appendix F.2 Multiple dilution likelihood for imperfect assay
sensitivity and specificity

Modifications to Web Appendix F.1 for different and multiple dilutions are
straightforward. First, consider assay data from a single dilution level of
$u\times 10^{6}$ cells per well, for some $u>0$. The imperfect assays
likelihood at this dilution level is proportional to
$\widetilde{L}(\boldsymbol{\tau}|\boldsymbol{W}^{*},\boldsymbol{Z}^{*},\boldsymbol{R},u)=L(u\boldsymbol{\tau}|\boldsymbol{W}^{*},\boldsymbol{Z}^{*},\boldsymbol{R})$,
following the same $\boldsymbol{\tau}=\boldsymbol{\lambda}/u$ substitution
from Section 2.6 in the main text. Similarly, the imperfect assays likelihood
for $D$ distinct dilution levels, $D\in\\{1,2,\dots\\}$, is proportional to

$\displaystyle\prod_{d=1}^{D}\widetilde{L}(\boldsymbol{\tau}|\boldsymbol{W}^{*(d)},\boldsymbol{Z}^{*(d)},\boldsymbol{R}^{(d)},u^{(d)}),$
(S.22)

where $\boldsymbol{W}^{*(d)}$, $\boldsymbol{Z}^{*(d)}$,
$\boldsymbol{R}^{(d)}$, and $u^{(d)}$ denote the vector of imperfect QVOA
results, matrix of imperfect UDSA results, vector of complete data indicators,
and dilution level, respectively, for the $d$th dilution. As in Section 2.6
the imperfect assays MLE for the vector of DVL-specific IUPMs is the value
$\widehat{\boldsymbol{\tau}}$ that maximizes (S.22), and the corresponding MLE
for the IUPM can be calculated as
$\widehat{T}=\sum_{i=1}^{n}\widehat{\tau}_{i}$.

###  Web Appendix F.3 Simulations for imperfect assay sensitivity and
specificity

The performance of the imperfect assays MLE was investigated in a simulation
study. Following the steps outlined in Section 3.1, true assay results $W$ and
$\boldsymbol{Z}$ were simulated for $M=12,24$, or 32 replicate wells at a
single dilution level. A dilution of $u=1$ million cells per well was chosen
with a true IUPM of $1$, such that $T=\Lambda=1$. The IUPM was allocated at a
constant rate across $n^{\prime}=6$ underlying DVLs, such that
$\lambda_{i}=1/6$, $i\in\\{1,\dots,6\\}$. For comparison, the “‘perfect
assays” IUPM estimator from the main text was included, which assumed 100%
sensitivity and specificity for both the QVOA and UDSA.

Three sets of simulations were conducted at a single dilution level. In the
first set, combinations of QVOA and UDSA sensitivities between $80\%$, $90\%$,
and $100\%$ were considered, with specificity fixed at $90\%$ for both assays.
In the second set, the assays’ specificities were varied while their
sensitivity was fixed at $90\%$ and $q=1$. In the first two sets of
simulations, all wells found to be HIV-positive using the imperfect QVOA would
undergo deep sequencing (i.e., $q=1$). In the third set, the proportion of
wells testing positive that underwent deep sequencing varied between
$q=0.5,0.75$, and $1$. In these simulations, the QVOA and UDSA had 90%
sensitivity and 90% specificity.

Figures 2 and 3 show the empirical distributions of the imperfect assays MLE
and the perfect assays MLE of the IUPM under varying assay sensitivities and
specificities, respectively. As expected, the perfect assays MLE was biased
when the perfect sensitivity and specificity assumption was violated. This
bias increased as the sensitivity and specificity decreased and persisted as
the number of replicate wells $M$ increased. Across all settings considered,
the imperfect assays MLE had very little bias, even as the amount of missing
deep sequencing data increased (Figure 4).

Figure 2: Empirical distributions of the imperfect assays MLE and the perfect
assays MLE of the IUPM at a single dilution level. QVOA and UDSA sensitivities
varied, but both assays had 90% specificity. Two replicates where the
imperfect assays MLE was $>10$ were excluded from the plot. Figure 3:
Empirical distributions of the imperfect assays MLE and the perfect assays MLE
of the IUPM at a single dilution level. QVOA and UDSA specificities varied,
but both assays had 90% sensitivity. One replicate where the imperfect assays
MLE was $>10$ was excluded from the plot. Figure 4: Empirical distributions of
the imperfect assays MLE and the perfect assays MLE of the IUPM at a single
dilution level. The proportion of imperfectly HIV-positive wells to undergo
deep sequencing varied, but both assays had 90% sensitivity and 90%
specificity. One replicate where the imperfect assays MLE was $>10$ was
excluded from the plot.

An additional set of simulations was conducted at multiple dilution levels.
Following the steps outlined in Section 3.2, true assay results at $D=3$
dilution levels were simulated with
$\boldsymbol{M}=(M_{1},M_{2},M_{3})=(6,12,18)$ replicate wells per dilution
level. The true IUPM of $1$ was constantly distributed across $n^{\prime}=6$
underlying DVLs, such that $\lambda_{i}=1/6$, $i\in\\{1,\dots,6\\}$. For the
imperfect assay results, combinations of QVOA and UDSA sensitivities between
$80\%$, $90\%$, and $100\%$ were considered, with specificity fixed at $90\%$
for both assays. The imperfect assays MLE continued to perform well for assay
data from multiple dilution levels, with very little bias across various error
settings (Figure 5).

Figure 5: Empirical distributions of the imperfect assays MLE and the perfect
assays MLE of the IUPM at multiple dilution levels. QVOA and UDSA
sensitivities varied, but both assays had 90% specificity. One replicate where
the imperfect assays MLE did not converge was excluded from the plot.

###  Web Appendix F.4 Considering imperfect assays in the HIV application

The HIV data from Section 4 of the main text were used to demonstrate how
imperfect assay sensitivity and specificity can be considered in practice.
Specifically, the IUPM estimates of Subject C13 were computed under 90%
specificity and varied sensitivity (the same settings considered for the
simulations in Figure 5) and compared to the IUPM estimates assuming perfect
sensitivity. Note the methods in Web Appendix F.3 above require per-well data
(i.e., $W_{j}^{*}$ and $Z_{ij}^{*}$). Unfortunately, the QVOA and UDSA results
were summarized by the lab so that such per-well level data was not available.
Therefore, well-level data consistent with the observed summary data were
simulated in order to illustrate application of the methods in Web Appendix
F.3 The resulting IUPM estimates and 95% confidence intervals for Subject C13
are shown in Figure 6. As as either the UDSA or QVOA sensitivity increased
from 80% to 100%, the IUPM estimates decreased and the corresponding
confidence intervals became narrower. This process could be repeated for all
subjects (i.e., source populations).

Figure 6: Estimated IUPM with 95% confidence intervals for Subject C13,
assuming imperfect assays. QVOA and UDSA sensitivities varied, but both assays
were assumed to have 90% specificity. The dashed line denotes the estimated
IUPM assuming perfect assays.