Real-time joint estimation of the size of the 2026 Bundibugyo virus outbreak in the Democratic Republic of the Congo
An outbreak of Ebola disease caused by Bundibugyo virus is ongoing in the Democratic Republic of the Congo, with exported cases in Uganda. Most infections are unreported, so the current size of the outbreak has to be inferred from the surveillance counts that are published. Published estimates use a subset of those counts under fixed delay, doubling-time and case-fatality assumptions. We instead fit every count the national situation reports publish in one posterior, refitted and released with each data update. One latent infection process by province, staged through explicit delays and stream-specific ascertainment, generates every stream. We evaluate each release by one-week-ahead forecasts scored against the data that arrive next and against a persistence baseline, by fits to single streams, and against the published estimates. Language-model agents drafted the model, its priors and the analysis under our direction, with each modelling decision taken by us. We read the development record to ask which events changed the estimates and forecasts and how each was found. At the release analysed the joint model estimates an outbreak larger than the confirmed count with the reproduction number near one, and its estimate sits below every estimate that the single streams imply on their own. Its forecasts beat persistence at most releases only for the onset curve, and its interval on size is no narrower than the priors on ascertainment and delays that the data do not move. In the record most model changes came from a person, most implementation defects were found and fixed by the agents, and the prospective evaluation caught defects that neither did. A joint estimate from published reports alone is feasible in real time and makes the agreement and conflict between surveillance streams explicit.
1 Centre for Mathematical Modelling of Infectious Diseases, London School of Hygiene & Tropical Medicine, London, United Kingdom
✉ Correspondence: Sam Abbott <sam.abbott@lshtm.ac.uk>
Introduction
Ebola disease caused by Bundibugyo virus has re-emerged in the Democratic Republic of the Congo (DRC), with exported cases in Uganda [1,2]. The number of infections is unknown while the outbreak is ongoing, and its size and growth drive the response, from bed and laboratory capacity to the risk of further exports [3]. Surveillance observes the reported fraction of infections after delays from infection to onset, confirmation and death, so the current size has to be inferred from the counts that are published.
The national situation reports carry many counts but no line list, and every count is report-dated and open to revision [1,4]. They report suspected and confirmed cases and deaths, laboratory throughput, isolation occupancy, treatment-centre flows, and an intermittently published curve of confirmed cases by onset date. Each stream observes the same epidemic through a different delay and a different ascertainment, so used alone each implies a different outbreak size.
Several estimates of the outbreak size have been published, but each uses a subset of these counts under fixed assumptions. McCabe et al. [5–7] estimated cumulative infections from the exported cases and the reported deaths under constant exponential growth, following the export-based approach of Imai et al. [8]. They report scenarios over fixed delay, doubling-time and case-fatality assumptions rather than one posterior. Chamla et al. [9] projected confirmed cases with a stochastic compartmental model fitted to the confirmed series alone. Genomic analyses bound the time of the most recent common ancestor and the early doubling time but not the size of the outbreak [10–12]. Earlier Bundibugyo outbreaks were small and give few outbreak-specific parameters, so the delay and severity assumptions every one of these estimates needs are drawn from other outbreaks or from judgement [13,14]. Renewal-based real-time tools such as EpiNow2 fit a case series with one secondary observation [15], and they are not designed for a set of report-dated aggregates that stop and change format.
In this paper we present a joint model of the outbreak size fitted to every count the situation reports publish, and its evaluation across four months of releases. The core idea is that one latent infection process, staged through explicit delays and stream-specific ascertainment, generates every published count, so all streams inform one posterior. We started from the approach of McCabe et al. [5], replaced it with a discrete renewal model in which the reproduction number varies over time [16,17], and then ran one renewal process per province. We release each update with the model definition, the draws and the forecasts. The model, its priors, its analysis and this text were drafted by language-model agents under our direction, with each modelling decision taken by us, and we treat that workflow as part of what we report. We read the development record for the events that changed the estimates or forecasts, how each was found and whether a person or an agent decided the response. We evaluate the model on the situation reports [1] by one-week-ahead forecasts against a persistence baseline, by fitting each stream alone, and against the estimates of McCabe et al. and Chamla et al. [5,9].
Methods
Data
The release analysed is results-2536, with data to 25 September 2026, and every current-model result below comes from its fits. DRC counts come from the daily situation reports of the Institut National de Santé Publique, except the 18 to 22 May suspected totals, which come from the WHO AFRO joint report [1,18]. Each value is read from the report PDFs by a language-model agent and re-read independently by a second agent, with disagreements settled against the source. The INRB-UMIE transcription is used only when the primary source is not publishing [4]. Exported cases enter as three dated detections and one dated death, frozen at the cut-off of the first situation-report model, and a fourth export reported afterwards is deliberately not fitted [2]. The curve of confirmed cases by symptom-onset date is digitised bar by bar from the figure the later reports carry. Genetic estimates of the time to the most recent common ancestor and of the early doubling time enter as priors, referred to below as the genetic bound [10,12]. Each stream is listed with its observation model, its first and last fitted report and the dates at which its reporting changed or was absent in Table 1.
The cumulative streams are fitted as the increments between successive reports, with a fitted level offset at each documented reclassification or harmonisation date. A stream that stops being published is frozen at its last report and stays in the likelihood.
Model development
The work began as a reproduction of the approach of McCabe et al. as one posterior rather than a set of scenarios, with exact delay convolutions in place of their closed-form approximations and a prior on every fixed quantity [5]. The national situation reports then became the data source, so the observation model moved to the per-report increments of running totals and added the confirmed series through a laboratory model that the throughput counts did not identify. A discrete daily renewal equation with a time-varying reproduction number replaced the latent process, chosen from five prototyped processes that included an explicit convolution, a stochastic renewal, a particle-filtered state space and a compartmental model [16,17]. Every remaining published stream was then added with its own observation model, so that each published count entered the likelihood on its own terms. One renewal process per province replaced the single national process when the reports began to carry counts by province that a national process could not use. The versions, the streams each of them fitted and the reason each gave way to the next are given in Table 2, and the estimate each released in Figure 3.
Current model
The full model definition and every prior are given in the Supporting Information (Methods). The model runs one renewal process for each of four provinces, Ituri, Nord-Kivu, Haut-Uele and one pooling the remaining affected provinces, coupled by a gravity kernel with same-day transfer, and every national stream is fitted against the sum over provinces. Each province \(p\) generates infections from its reproduction number \(R_{p,t}\) and the generation interval \(g\), with the log reproduction number the national trend plus a province deviation that sums to zero across provinces
\[ G_{p,t} = R_{p,t} \sum_{s \ge 1} I_{p,t-s}\, g_s, \qquad \log R_{p,t} = \log R^{\text{trend}}_t + \delta_{p,t}, \qquad \sum_p \delta_{p,t} = 0. \tag{1}\]
The national trend follows a weekly random walk, and the deviations revert toward zero at each knot on an orthonormal sum-to-zero basis and are scored through the per-province compositions of confirmed cases and deaths [19,20]. A single import grows through a cryptic phase to the renewal start, with the genetic bound flooring the outbreak age and the number of cryptic generations given a prior. Importation relocates a share of what each province generates, with \(K\) the gravity kernel and \(\varepsilon_{q,t}\) the share leaving origin \(q\), and national infections \(I_t\) are the sum over provinces
\[ I_{p,t} = \Bigl(1 - \varepsilon_{p,t} \sum_{q \ne p} K_{q,p}\Bigr) G_{p,t} + \sum_{q \ne p} \varepsilon_{q,t} K_{p,q}\, G_{q,t}, \qquad I_t = \sum_p I_{p,t}. \tag{2}\]
The generation interval \(g\) is a gamma distribution whose mean and standard deviation carry priors centred on the serial interval of Ebola virus disease in West Africa, with widths from a Bundibugyo cohort [21,22].
National infections convolved with the incubation period give daily onsets \(o_t\), and each stream convolves those onsets with its own sampled and double-censored delay, scales them by its ascertainment and adds any stream-specific term [23]. For the suspected cases the delay is the onset-to-report distribution \(f_{\text{rep}}\), the ascertainment is \(p_{\text{DRC}}\) and the added term is a non-BVD background rate \(\lambda_{\text{bg},t}\) on a weekly random walk. The increment between the reports on days \(d_{i-1}\) and \(d_i\) is scored with a negative binomial whose dispersion \(k_s\) is partially pooled across the four passive-surveillance streams
\[ c_t = p_{\text{DRC}} \sum_{s \ge 0} o_{t-s}\, f_{\text{rep},s} + \lambda_{\text{bg},t}, \qquad Y_i - Y_{i-1} \sim \mathrm{NegBinomial}\Bigl(\sum_{t = d_{i-1}+1}^{d_i} c_t,\ k_s\Bigr). \tag{3}\]
The model estimates infections and deaths to date, the reproduction number by province, the national growth rate, the case-fatality ratio, and the ascertainment of each surveillance system.
Forecast evaluation
Every release forecasts each still-reported stream one to four weeks ahead, and we score the one-week horizon once the next reports arrive. Forecasts are drawn from the fitted model by running it past the cut-off with the fitted parameters fixed, so each future count comes from its stream’s own likelihood, and the national forecast is split across provinces by the fitted compositions. The scored releases predate this method and projected each stream from its rate at the cut-off under the continued reproduction-number walk. We score with the continuous ranked probability score, where lower is better, on the count scale and on a log scale that stops the largest counts dominating [24,25]. We also report coverage of the 50% and 90% predictive intervals [26].
The persistence baseline carries the increment over the preceding window forward, with a spread drawn from the stream’s own past differences. It was checked to read only the reports available on the day the forecast was made, since the choice of baseline sets how much skill a forecast appears to have [27]. Streams no longer reported are excluded because their truth is zero by construction, which favours neither the model nor the baseline. Fits fixed at each release date are scored on one scale across reporting breaks, with the harmonisation step subtracted from the observed counts and the windows holding a reclassification dropped.
Comparison with single-stream fits and published estimates
The joint estimate is set against the same model fitted to one stream at a time with a single likelihood term, against the estimate of every earlier release, and against the published estimates. The estimate of each earlier release is read from that release’s published outputs. McCabe et al. are matched at their cut-offs by fits fixed at those dates, and Chamla et al. by a projection from a fit fixed at their calibration date [5,7,9].
Convergence
Each fit reports R-hat, bulk effective sample size and the number of divergent transitions for every parameter [28,29]. We sample with the No-U-Turn sampler (NUTS), with 2 chains of 1,000 draws after 500 adaptation steps, a maximum tree depth of 10 and a target acceptance probability of 0.80 [30]. Each chain starts from the first of a batch of prior draws whose log density is at or above the batch median, since a start in the prior tail can trap a chain. The build fails if the joint fit exceeds an R-hat of 1.1 or a divergent fraction of 5%, or falls below a bulk or tail effective sample size of 25, and it warns at a tighter tier. Builds after the release analysed also fit the joint model to a dataset simulated from one of its prior draws, with the same sampler settings, and check that the posterior recovers the generating values and the simulated future.
Development process
Language-model agents drafted the code, the priors and the analysis, and we directed the work from issues and pull-request reviews and took each modelling decision. We worked this way because one modeller’s time was set against a live outbreak, because we aimed to fit every published count, and to learn what the workflow could and could not do. An agent opened a pull request for each change, and the lead author reviewed it locally before anything was recorded on the request. Later an automated reviewer commented on each pull request, and the author accepted or rejected each of its findings on the request. A change that did not converge or did not fit its stream was not merged, judged by review of the diagnostics and later by a check in continuous integration. The first releases were signed off by every author on a release issue opened for the purpose. Between changes the agents ran wide searches over parameterisations, samplers and gradient backends, and most of these measured as null and were not merged.
Beyond the review of single changes, one check applied to the transcribed data and another to the released model as a whole. The second read of each transcribed value, described under Data, checked the transcriptions, and by the lead author’s rule every fitted stream advanced on every data update so that no stream was left behind. The prospective evaluation, which scored each release against the data that arrived next, checked the released model on data that no review or test had seen.
Analysis of the development record
We read the development record for the events that changed the released estimate, the reproduction-number trajectory or the forecasts, and for the defects in released outputs, to establish how each was found and who decided the response. The sources were the issues, the pull-request bodies and review comments, the automated reviewer’s comments and the release notes, with the account behind each item recorded. An item opened under a human login counted as agent-written where its body carried the marker the agents leave, since the human accounts also ran agent sessions. We classified each event by kind, as a model change, a data change, a data-quality decision or a defect, and by the route through which it was first recorded. We also recorded who decided the response, a person, an agent or the two jointly, and kept the evidence link for every row (Table 3). An agent counted as the finder when the first record was an issue, pull-request body or comment from an agent account, and as the decider when the record held only agent text and at most an approval without a body. Where releases fell on both sides of an event we recorded the two released medians, marked whether other changes shared that window, and estimated no effect ourselves.
From the same record we counted commits and merged pull requests per week by account, human-authored issues, review submissions and inline comments, and the automated reviewer’s submissions and inline comments (Figure 3 B). The record cannot show decisions taken off GitHub, and the lead author’s local first-pass reviews appear on it only as an approval.
Implementation
The model is written in Julia with Turing and differentiated with Mooncake [31,32]. The daily convolution, renewal and observation kernels carry hand-written reverse-mode rules, each checked against central finite differences and against the backend’s own derivation of an unregistered copy. The agents wrote those rules and a precompile workload that compiles the fit the report runs, and the benchmarks recorded for both are in the Supporting Information (Gradient and build benchmarks). Each data update refits every affected model in continuous integration, renders the report and publishes a versioned release with the model definition, thinned draws and forecasts.
| Stream | Observes | Likelihood | First | Last | Breaks and gaps |
|---|---|---|---|---|---|
| Exported cases | Cases detected in Uganda, by detection date | Poisson, delayed by onset to detection abroad | 11 May | 23 May (frozen) | none |
| Deaths among exports | Deaths among exported cases | Poisson on the dated death, delayed by onset to death | 14 May | 14 May (frozen) | none |
| Suspected cases | Cumulative suspected cases, per report | Negative binomial increments with a non-BVD background | 18 May | 26 May (frozen) | none |
| Suspected deaths | Cumulative suspected deaths, per report | Negative binomial increments with a background | 18 May | 26 May (frozen) | none |
| Daily new suspects | Suspects notified in the previous 24 h | Negative binomial | 4 June | cut-off | 6 August absent |
| Daily suspected deaths | Suspected deaths in the previous 24 h | Negative binomial | 7 June | 11 July (frozen) | none |
| Confirmed cases | Cumulative laboratory-confirmed cases, per report | Beta-binomial of the specimens analysed where a denominator is published, negative binomial increments otherwise | 14 May | cut-off | 22 July |
| Confirmed deaths | Cumulative confirmed deaths, per report | Negative binomial increments through the death-side pipeline | 14 May | cut-off | 22 July |
| Laboratory total | Cumulative samples analysed | Negative binomial | 23 May | 28 May (frozen) | none |
| 24 h analysed volume | Samples analysed in the previous 24 h | Negative binomial, anchors positivity | 1 June | cut-off | none |
| Isolation occupancy | Patients in isolation on the report day | Censored negative binomial against a supply-limited bed model | 1 June | cut-off | 9 June, 19 June, 14 July, 11 August |
| Bed capacity | Isolation beds available | Negative binomial of the implied bed count around a log-scale random walk on capacity | 9 June | cut-off | none |
| Recovered | Cumulative confirmed recoveries | Negative binomial increments | 6 June | cut-off | none |
| Treatment-centre flows | Admissions, in-care deaths, rule-outs and abscondings | Negative binomial increments with an in-care case-fatality ratio | 13 June | 2 August (frozen) | none |
| Onset curve | Confirmed cases by symptom-onset date, per report | Student-t reporting triangle with a fitted read error per bar | 12 July | cut-off | absent where a report was not published or carried no onset figure |
| Province compositions | Confirmed cases and deaths by province, with samples analysed per head as an ascertainment covariate | Beta-binomial stick-breaking split of the national count | 6 June | cut-off | none |
| Version | Released | Latent process | Streams | Reason for the next version |
|---|---|---|---|---|
| v1.0.0 to v1.2.0 | 20 to 27 May | Continuous-time constant exponential growth with exact delay convolutions | Exports, deaths among exports, suspected cases and deaths as cumulative totals | The situation reports gave per-report increments and the suspected totals were withdrawn |
| v1.3.0 | 9 June | As above | Per-report increments; confirmed cases and deaths through a laboratory queue | The laboratory model was not identified by the throughput counts and the joint estimate sat above every single stream |
| v1.4.0 to v1.18.0 | 9 June to 9 September | Daily renewal with a weekly random walk on the reproduction number and a cryptic seeding phase | Every stream in Table 1 except the province compositions, added over the window | To fit the counts by province that the reports had begun to carry |
| v2.0.0 to v2.1.0 | 16 to 18 September | One renewal process per province, coupled by importation | Adds the province compositions | The first provincial fit did not converge, and centring the background and bed-capacity walks restored it |
| v2.2.0 | Unreleased | As above, with the province deviations on a sum-to-zero basis and hand-written gradient rules | As above | The version analysed here; main has since added susceptible depletion, province splits of the laboratory and isolation counts, a prior on the size of the province deviations and a narrower generation-interval prior |
Results
The first release
At the first release the joint posterior of cumulative infections sat above both headline scenarios of McCabe et al. [5], a median of 972 (90% interval 478 to 2,437) against 313 and 501 (Figure 1 A). The single-stream fits overlapped one another and the joint fit, with the suspected-death fit the narrowest and the fit to deaths among exports the widest (Figure 1 B). The latent curve placed most of the seeding density in the months before the first report, and both fitted DRC counts sat below the median curve at the cut-off (Figure 1 C).
Outbreak size and dynamics at the release analysed
At the data cut-off of the release analysed the joint model estimated 16,000 cumulative infections (90% interval 11,000 to 23,000), 2.0 times the 7,989 confirmed cases reported by then. The posterior daily onsets rose over the first two months of reporting, levelled off from July, and eased towards the cut-off (Figure 2 A). The national reproduction number at the cut-off was 0.92 (90% interval 0.67 to 1.2), with 27% of draws above one (Figure 2 B). The daily growth rate at the cut-off was -0.03 to 0.02 (90% interval), the case-fatality ratio was 50% (90% interval 37% to 61%) and the first infection was placed 182 to 209 days before the cut-off (90% interval). Every provincial reproduction number straddled one at the cut-off, and the probability that it was above one was 24% in Ituri, 32% in Nord-Kivu, 44% in Haut-Uele and 74% in the remaining provinces (Figure 2 B).
The posterior on suspected-case ascertainment moved below its prior median of 0.75 to 0.58 (90% interval 0.37 to 0.82), and export ascertainment moved from 0.75 to 0.64 (90% interval 0.39 to 0.88) (Figure 2 C). The 90% interval of the fraction of suspects tested, 0.73 to 0.97, sat almost entirely above its prior median of 0.74, and the onset-report ascertainment was 0.28 to 0.59 (90% interval) (Figure 2 C).
The fit that produced these estimates passed the convergence gate, with a maximum R-hat of 1.032, a minimum bulk effective sample size of 73 and a divergent-transition count of 0 over 2 chains of 1,000 draws. The forecast made at the previous release gave 662 new confirmed cases (90% interval 485 to 875) and 291 new confirmed deaths (90% interval 226 to 375) for the week to the cut-off (Figure 2 D). The reports carried 375 cases and 176 deaths over that week, so both streams fell below the 90% interval of the forecast (Figure 2 D).
The release before the one analysed did not converge and was published with a convergence warning. Convergence was restored at the next release, and the released median stayed inside the interval of the release that failed (Figure 3 C).
Estimate evolution across releases
The released estimate rose by more than an order of magnitude, from 972 at the first release to 13,314 at the last tagged release (Figure 3 C). The largest step came with the move to situation-report increments and the laboratory model, where the median rose by a factor of 3.2 between the last two releases of the closed-form model. Over the same window the number of fitted streams rose from four to eleven and most commits came from the agent account (Figure 3 A and B). Under the provincial model the two released estimates plateaued, and the release analysed, which postdates the last version tag and is not drawn, sat above them.
Development record
| Date | Event | Effect | Found by | Decided by | Releases live |
|---|---|---|---|---|---|
| 26 May | Data source moved to the INSP situation reports | Released median 925 to 1,391 (v1.1.0 to v1.2.0) | human review | human | |
| 6 June | Joint laboratory-queue fit sat above every single-stream fit | Released with a warning banner, then replaced by the renewal model | agent | human | 1, from v1.3.0, fixed in v1.4.0 |
| 9 June | Renewal model with a weekly random walk on R replaced the constant-growth model | Released median 4,509 to 4,161 (v1.3.0 to v1.4.0) | human review | human | |
| 10 June | 24h analysed laboratory volume fitted to anchor positivity | Released median 4,161 to 4,490 (v1.4.0 to v1.5.0) | human data check | human | |
| 10 July | Molecular-clock and growth priors updated to the outbreak-specific genomic analysis | Released median 4,380 to 4,791 (v1.8.0 to v1.9.0) | human data check | human | |
| 4 August | Growth-rate prior widened | Released median 7,446 to 8,672 (v1.12.0 to V1.13.1) | human data check | human | |
| 10 August | Brief-format reports dropped the treatment flows and the daily suspects, both frozen | Released median 8,672 to 9,147 (V1.13.1 to v1.13.2) | agent | human | |
| 26 August | Five defects in the one-week-ahead forecast, among them the last innovation compounded over the horizon | 90% forecast coverage sat at one until the fix | prospective evaluation | joint | 13, from v1.4.0, fixed in v1.15.0 |
| 28 August | Reporting-break corrections applied in the fit but not in the scoring | Frozen-fit forecasts scored an order of magnitude out | prospective evaluation | human | 8, from v1.11.0, fixed in v1.17.0 |
| 3 September | Scoring fix carried confirmed but not occupancy break days onto frozen fits | Bed-stream scoring across breaks | automated reviewer | agent | 8, from v1.11.0, fixed in v1.17.0 |
| 3 September | Bed shortfall could go negative under a positive offset | Forecast bed-shortfall column | automated reviewer | agent | 0, found before release |
| 4 September | Onset triangle given a per-scan level error | Released median 12,910 to 12,100 (v1.16.0 to v1.17.0) | agent | joint | |
| 14 September | Analysed laboratory volume capped below the modelled suspect inflow, against the published methods | Released fits before the provincial model carried the cap | agent | joint | 16, from v1.5.0, fixed in V2.0.0 |
| 15 September | Four-province model became the headline | Released median 13,771 to 14,578 (v1.18.0 to V2.0.0) | human review | human | |
| 17 September | Provincial release published with an unconverged headline fit | Released median 14,578 to 13,314 (V2.0.0 to v2.1.0) | agent | agent | 1, from V2.0.0, fixed in v2.1.0 |
| 17 September | Matched forecast scoring keyed without the made date | Frozen confirmed-death skill 3.18 to 1.76 | agent | agent | 11, from v1.11.0, fixed in v2.1.0 |
| 18 September | Provincial confirmed composition missed the report leg of the delay | Provincial split attributed to earlier days than the national total | agent | agent | 1, from V2.0.0, fixed in v2.1.0 |
| 21 September | Sensitivity page still described the headline as three provinces | Released report text | agent | agent | 2, from V2.0.0, fixed in v2.2.0 |
| 24 September | Onset digitiser read 10 to 23% short of the printed total and mis-dated bars | Onset curve data from v1.16.0 to v2.1.0 | agent | joint | 5, from v1.16.0, fixed in v2.2.0 |
| 24 September | Comparison table printed the current confirmed total as observed at the calibration date of Chamla et al. | Released sensitivity page, fix not yet released | agent | human | 17, from v1.7.0, fix unreleased |
The record held 49 events, of which 17 were model changes, 16 defects, 11 data-quality decisions and 5 data changes (Figure 3 D). The agents first recorded 26 of them and a person 18, with the prospective evaluation 2, the automated reviewer 2 and no route recorded for the remaining 1. The routes differed by kind, with 11 of the 17 model changes first raised by a person and 10 of the 16 defects first reported by the agents. A person decided 23 of the events, an agent 17 and the two jointly 9, and the agent count is an upper bound because the lead author’s local first-pass reviews appear on the record only as an approval (Figure 3 D).
The largest moves in the released median each straddled several changes merged between the same two releases, so none can be attributed to one event (Figure 3 C). The largest, reported above, followed the reclassification of the suspected counts, which an agent reported and the lead author answered by freezing those streams and fitting the confirmed ones alongside the laboratory model. The next, a rise by a factor of 2.9, came at the release that gave a reclassification of the isolation census a fitted level offset, proposed by an agent and approved without comment. The third, a fall by a factor of 2.0, came at the release in which an agent widened the prior on the non-BVD background with no human comment on the record, alongside the isolation stream and an earlier start to the random walk.
The events the agents found themselves included the reclassification of the suspected counts, the lag in the export stream, the streams lost when the reports moved to a brief format, the unconverged provincial release and a chain frozen at its start by a fixed seed. Three agent explanations were later withdrawn or overturned, a diagnosis on which one onset block was excluded and then restored, a fetch failure called permanent and then transient, and a factor pinning the provincial reproduction numbers to the national trend that was defended with a prior simulation and then removed (Table 3).
The 2 findings of the automated reviewer both came in the pull request that corrected the scoring, a break-day correction that had not reached the bed stream and a shortfall that could go negative, and the agent accepted the first and declined the second (Table 3). The reviewer declined 3 further pull requests over its line cap. The 2 defects found only by the prospective evaluation were the compounded forecast innovation and the scoring without break corrections, both read by the lead author from the forecast panels of the validation page (Table 3). The defects that reached the forecasts were carried by as many as 13 releases before the scoring exposed them (Table 3).
Over the window the agent accounts made 542 commits and had 422 pull requests merged, against 59 commits and 21 merged requests from the human accounts (Figure 3 B). The human accounts submitted 556 reviews with 183 inline comments and the automated reviewer 175 reviews with 244 inline comments. Over the 19 weeks of the record that is 29 reviews and 9.6 inline comments a week from the human accounts against 28 agent pull requests, and it is the only measure of a person’s time the record holds. The window carried 153 releases, 8.1 a week, each refitting every affected model, against one recorded joint fit wall-clock of 320 minutes.
Forecast evaluation
The joint one-week-ahead forecast scored worse than persistence on most stream-release pairs, beating the baseline on 11 of 69 pairs overall and on 9 of 53 at the live releases (Figure 4 A). By stream it beat persistence on 1 of 20 releases for confirmed cases, 0 of 20 for confirmed deaths, 3 of 18 for isolation beds, 6 of 8 for onset reports and 1 of 3 for recovered. The 90% interval of the joint forecast covered 67 of 69 observations (Figure 4 A). Every scored forecast predates the correction of the forecast defects reported in the Discussion, so the intervals scored here are wider than the model now produces.
Comparison with single-stream fits and published estimates
The single-stream fits disagreed on cumulative infections at the cut-off by a factor of 13.16, from 17,437 for the isolation beds fit to 229,543 for the exports fit (Figure 4 B). The joint estimate of 15,795 sat below the median of every single-stream fit (7 of 7), and its 90% interval of 11,000 to 23,000 overlapped only the isolation beds, onset reports and exports intervals (Figure 4 B).
At the second cut-off of McCabe et al. [5,7] the frozen joint interval overlapped the upper half of their back-calculation scenarios and the upper parts of their geographic-spread intervals (Figure 4 C). At their third cut-off the frozen joint interval overlapped both sets of scenarios, and the back-calculation scenarios sat in its lower half (Figure 4 C).
The joint fit frozen at the calibration date of Chamla et al. [9] projected 1,474 confirmed cases (90% interval 1,022 to 2,494) by their week-12 date, against their central 990 (709 to 1,293) (Figure 4 D). The observed count of 1,155 lay above their central projection and inside both their interval and the joint’s (Figure 4 D).
Discussion
In this paper we have presented a joint model of the size of the 2026 Bundibugyo virus outbreak fitted to every count the national situation reports publish. At the release analysed the joint estimate sat below every single-stream fit. The one-week-ahead forecasts rarely beat a persistence baseline, and the onset reports were the only stream on which they beat it at most releases.
The approach has several strengths. First, one posterior over every published stream shows where the streams agree and where they disagree, and resolves that disagreement inside the model. Under the renewal model that disagreement appears to be absorbed by the ascertainment and background terms, while under the laboratory model it inflated the size estimate above every single stream. Second, fitting the cumulative counts as increments with a level offset at each documented change handles revisions and format changes rather than ignoring them. Finally, the whole model definition lives in one versioned document that others can review and adapt, and it is refit and released with each data update.
The main limitation is that the delays, the case-fatality ratio and the laboratory parameters rest on priors from other outbreaks or from judgement. The data do little to move them, so the size interval inherits their width and any error in them moves the estimate the same way. Almost every count is report-dated, so the timing of the epidemic is recovered through those assumed delays rather than observed. The release analysed carries the corrected generation-interval prior, but the fitted mean still sits at its lower edge, which feeds into every reproduction number reported here. The reproduction of the first model has not been checked against the original authors’ code. The data are aggregate national and provincial counts, with no line list in the released model.
The streams share one case pool but are treated as conditionally independent given the latent infections, which may understate the uncertainty in the joint estimate. The fitted streams may look sub-exponential while the reproduction number sits above one, so the model may misstate growth over the window.
Each published estimate answers a different question from a different subset of the published counts. The reproduction of the method of McCabe et al. [5] at their own fixed inputs recovers their headline estimate, and the joint posterior sits above both of their headline scenarios and overlaps only the upper half of their back-calculation scenarios at the matched cut-off. What the joint model adds is the streams that inform ascertainment, which their fixed assumptions leave outside the posterior. It adds to the projection of Chamla et al. [9] a forecast drawn from every stream rather than from the confirmed series alone. Genomic estimates bound the outbreak age and the early growth rate, and they enter the joint model as priors that the surveillance data then update [12]. More generally, the model belongs to the class of multi-stream renewal models with nowcasting components, applied here to report-dated aggregates rather than to line lists [15,33].
Changes in reporting, rather than changes in the epidemic itself, were the main source of change to the model. Most were absorbed by a fitted level offset or by freezing the affected stream, and the latent process was replaced only when a stream could not be fitted at all.
The record shows which decisions fell to a person and which to the agents, and the split follows the kind of event (Figure 3 D). Most model changes were raised and decided by a person, most implementation defects were first reported and fixed by the agents, and the defects that neither found were caught by the prospective evaluation (Table 3). The data-quality decisions, where a stream was frozen, given a level offset or excluded, went to a person in most cases, and in one a person overruled the agent’s reading of the series. In this project a person took each change of model structure, data rule and prior and read the diagnostics, while the agents carried the implementation. The person holding those decisions here had to read sampler diagnostics, weigh a prior against the published literature and know the surveillance system well enough to see a reclassification in a series. Working this way also made engineering feasible that we would not have attempted alone, among it the hand-written adjoints, the wide performance searches and the release tooling. We did not log token use or wall-clock time, so we cannot price this way of working here.
Several defects were of a kind a human author would be unlikely to write, and each passed the tests and the in-sample checks of the release that carried it. The forecast extrapolated the last random-walk innovation as a fixed slope, so its intervals widened with the horizon and covered every observation. Reporting-break corrections were applied in the fit but not in the scoring, so forecasts from fits fixed at each release date scored an order of magnitude out. The posterior predictive was generated without the provincial arguments, so every case-driven stream came out short while the background-driven streams were unaffected. A chain initialised from a fixed seed sat on a prior tail and returned the same R-hat on successive datasets, and a generation-interval prior documented as matching a published interval was much wider than that interval. The last two were caught in review of the pull request that carried them, and the record does not say how the posterior predictive defect was found.
Ranked by what they caught here, the prospective scoring came first, because it found defects that no review, test or automated reviewer had seen, and review of the pull request came next. The automated reviewer left more inline comments than the human accounts and had both its findings in one pull request, so it was marginal here, and the convergence gate in continuous integration postdates the one released fit it would have stopped. We would not delegate the choice of a prior again, because an agent widened the prior on the non-BVD background with no comment on the record and moved the released median. We would also not delegate the structure of the latent process, or an agent’s account of why a fit behaves as it does, since a factor pinning the provincial reproduction numbers was defended with a prior simulation and then removed in review. Whether the same division would hold for another team, another outbreak or a later stage of this one, the record cannot say.
A health-zone model conditioned on the provincial posterior is in development. Outbreak-specific delays could be estimated from a line list, following Akilimali et al. [22]. Revision could be treated as an observation process rather than as a break. A retrospective evaluation against the final counts could follow once the outbreak has ended.
The model has kept expanding and the experiments that shaped it were ad hoc, so future work could run formal experiments across settings and outbreaks. Our own infectious disease modelling workflow [34] could be applied more deliberately to test whether it improves the efficiency of agent work, measured as token use per accepted change. How agent-driven modelling should be directed is open, and on sharing the work this project can say only that one lead author held every modelling decision while the co-author reviewed the releases and this paper.
A joint estimate of outbreak size can be produced in real time from published situation reports alone. Its value lies in making the agreement and the conflict between surveillance streams explicit rather than in any one number. Its interval separates an outbreak twice the confirmed count from one ten times it, and the forecasts covered nearly every observation, so it supports the scale of a response and not its week-to-week planning. Agent-drafted modelling of this kind can work at outbreak speed when a person holds the modelling decisions. We aim to support real-time size estimation from published surveillance in future outbreaks.
Funding
This study is funded by the National Institute for Health and Care Research (NIHR) Health Protection Research Unit in Health Analytics & Modelling, a partnership between the UK Health Security Agency, Imperial College London and the London School of Hygiene & Tropical Medicine (grant code NIHR207404). The views expressed are those of the authors and not necessarily those of the NIHR, UK Health Security Agency or the Department of Health and Social Care.
Acknowledgements
We thank Katharine Sherratt and Samuel P. C. Brand for their contributions to the early versions of this work and their review of it, Neil Ferguson for suggesting the genetic bound on the seeding time, the Institut National de Santé Publique for publishing the situation reports and the INRB-UMIE team for their transcription, and the contributors who raised issues on the repository.
Data and software availability
AI disclosure
Claude Code sessions (https://www.anthropic.com/claude) drafted the model code, priors, tests, analysis text and this manuscript, transcribed the situation reports, and reviewed pull requests through an automated reviewer. The commit trailers record Claude Sonnet 4.6, Opus 4.8, Fable 5, Opus 5, Sonnet 5, Fable 5.1 and Opus 5.5, and the model behind the automated reviewer is not recorded. The authors set the modelling decisions, reviewed and revised all of it, and take responsibility for the content. Every quoted number is read from a released output or from the pull request that recorded it, rather than typed.