MATERIALS AND METHODS

Ethical approval

The surgical procedures (spaying and neutering) were performed by a team of qualified veterinarians registered with the Sri Lankan Veterinary Council. Since only routine procedures were performed and no experimental interventions involving live animals were undertaken, formal ethical approval for the surgeries was not required. Owned animals were sterilized only after obtaining the owners’ consent. Free-roaming animals were sterilized with the approval and supervision of the relevant local government authorities, namely the Anuradhapura Municipal Council and the Manupa, Nanupa, and Mihinthale Pradeshiya Sabhas. Free-roaming animals were captured, vaccinated, sterilized, and returned to their original habitats after adequate recovery in accordance with internationally recognized CNVR protocols [37] and the field clinic guidelines of the Department of Animal Production and Health. Animals that were sick, pregnant, lactating, very young, or very old were excluded from sterilization.

The questionnaire survey, which was conducted to determine the proportion of pets that were sterilized and vaccinated against rabies, was preceded by a detailed verbal explanation to pet owners regarding the importance of the survey, and written informed consent was obtained through signed response forms. Prior written permission to conduct the survey was obtained from the Provincial Commissioner of North Central Province. The survey was conducted in accordance with the principles of the Declaration of Helsinki.

Study period and location

The questionnaire survey was conducted from September to October 2025. The dog population control model program was initiated in May 2020 at the Municipal Veterinary Office of Anuradhapura, with the initial activities concentrated in the ancient city, which is characterized by a high density of free-roaming dogs. Based on observations made during 2020, the CNVR program adopted a more systematic approach beginning in March 2021. The Anuradhapura Municipality was divided into three circular zones according to demographic characteristics (Figure 1).

The first zone (8°20′26.01″N–8°18′29.03″N and 80°24′17.55″E–80°25′21.81″E) had a diameter of 2.6–2.9 km and an area of 5.94 km², encompassing the commercial city, consisting mainly of shops, offices, and transport hubs. The dog population in this region consisted predominantly of strays (free-roaming dogs without owners or responsible caregivers), representing the less catchable fraction.

The second zone, with a diameter of 4.3 km and an area of 8.58 km², was located between 8°21′10.68″N and 8°18′45.37″N and 80°23′23.64″E and 80°25′59.79″E and surrounded Zone 1. This zone primarily covered the ancient city and included only a limited number of residential areas, where dumping of unwanted puppies was frequently observed.

The third zone, covering 18.66 km² with a horizontal diameter of 6 km and a vertical diameter of 7 km, extended between 8°22′54.77″N and 8°18′9.05″N and 80°23′7.59″E and 80°26′21.71″E. This zone surrounded the first two zones and consisted mainly of residential areas. The dog population comprised owned and community dogs (dogs without a single owner but collectively cared for by community members), the majority of which were free-roaming. The boundary of this zone, represented by a red dashed line, approximated the municipal boundary of Anuradhapura. The total area of the municipality is 36.32 km², and the three circular zones collectively covered 33.18 km², corresponding to 91.35% coverage. The peripheral areas excluded from the third circle were included in the buffer zone.

Owners brought their dogs to the clinics for sterilization, whereas teams of dog catchers captured free-roaming stray and community dogs for ARV and sterilization. Both male and female dogs that had reached puberty were selected for sterilization without intentional sex bias. In contrast, owners tended to prioritize the sterilization of female animals.

The geographic coordinates of sterilization centers and catching locations were recorded using Google Earth Pro (±5 m accuracy), and all maps were developed using QGIS version 3.44.5.

During 2021 and 2022, sterilization centers and catching activities were confined to the three zones within the Anuradhapura Municipal Council boundaries. However, no marked reduction in the number of unsterilized free-roaming dogs or the dumping of new litters was observed by the end of 2022, highlighting the necessity of "sealing the boundary to cut down the supply." Consequently, the project boundary was extended by 3 km beyond the municipal boundary, resulting in a circle with a diameter of 10 km and an area of 45.36 km², excluding the first three zones. This area was designated as the buffer zone, represented by a black dashed line in Figure 1. The buffer zone was prioritized during 2024 and 2025 to minimize the influx of intact adult dogs and the dumping of puppies into Zones 2 and 3.

Figure 1: Distribution of sterilization centers and demarcation of zonal boundaries within the Anuradhapura Municipality from 2021 to 2025.

The zonal demarcation relied heavily on anthropogenic features such as commercial districts and residential areas, which are known to influence dog behavior [38]. Furthermore, municipal boundaries had to be considered because of state supervision. The actual land area of each zone was smaller than the calculated values due to the presence of large water reservoirs. The 10-km-diameter buffer zone (3 km beyond the municipal boundary) was selected because its anthropogenic characteristics and human population density are similar to those of Zone 3. Areas beyond 3 km are more sparsely populated and have a lower density of free-roaming dogs. Moreover, the roaming distance of Sri Lankan dogs rarely exceeds 2 km, making a 3-km buffer sufficient to minimize immigration into the city. Edge effects were not considered in the calculations.

The Anuradhapura Municipality is surrounded by three Pradeshiya Sabhas (local government areas), namely Nuwaragam Palatha Central (Manupa), Nuwaragam Palatha East (Nanupa), and Mihinthale. Regions within and beyond the buffer zone belong to one of these administrative divisions. Project activities were extended to these neighboring areas by 2022 as an additional precaution against the influx of dogs into the city.

Figure 2 illustrates the distribution of sterilization centers within the boundaries of the three Pradeshiya Sabhas and the Anuradhapura Municipality. Manupa is shown in yellow, Nanupa in cyan, and Mihinthale in pale green.

Table 1A summarizes the time frame and geographical coverage of the sterilization programs. Table 1B presents the numbers of anti-rabies vaccinations and sterilizations performed by the Municipal Veterinary Office from 2014 to 2019 for comparison.

The human population of the Anuradhapura Municipality is currently 63,276, and based on the national dog-to-human ratio of 1:8, the estimated dog population is approximately 7,910.

Conducting at least four sterilization programs annually, with intervals not exceeding 3 months, considerably reduced the probability of females that had escaped sterilization during previous campaigns becoming pregnant, because the reproductive cycle of dogs in Sri Lanka typically lasts approximately 6 months.

Outreach strategy and catching process

Areas with a high density of stray and free-roaming dogs within the municipality were identified during preliminary surveys conducted in 2020. Collaboration with voluntary dog feeders facilitated the identification of locations requiring special intervention.

Figure 2: Distribution of the sterilization centers within the Manupa, Nanupa, and Mihinthale Pradeshiya Sabha boundaries.

Year | The number of mass sterilization programs conducted per year | Total number of catching days per year | Zones covered | Anti rabies vaccination for dogs | Number of sterilizations
2020* | 2 | 8 | Preliminary events before the demarcation of the zones* | 777 | Male dogs: 94
Female dogs: 575
All dogs: 669
2021 | 6 | 54 | 1,2,3 | 5742 | Male dogs: 927
Female dogs: 1789
All dogs: 2716
2022 | 6 | 47 | 1,2,3, Buffer, Manupa, Nanupa, Mihinthale | 8128 | Male dogs: 765
Female dogs: 1862
All dogs: 2627
2023 | 3 | 40 | 1,2,3, Buffer, Manupa, Nanupa, Mihinthale | 6457 | Male dogs: 487
Female dogs: 1320
All dogs: 1807
2024 | 4 | 39 | 1,2,3, Buffer, Manupa, Nanupa, Mihinthale | 5334 | Male dogs: 266
Female dogs: 1110
All dogs:1376
2025 | 3 | 36 | 1,2,3, Buffer, Manupa, Nanupa, Mihinthale | 5528 | Male dogs: 248
Female dogs: 891
All dogs: 1139
Total |  |  |  | Male dogs: 2787
 |  |  |  | Female dogs: 7547
 |  |  |  | All dogs: 10,334

Year | ARV | Sterilizations
2014 | 3609 | 1467
2015 | 3673 | 1810
2016 | 4578 | 1159
2017 | 2924 | 1401
2018 | 2654 | 1509
2019 | 3206 | 1633

Two dog-catching teams, each comprising three members, were assigned two predetermined routes for each zone throughout the study period, regardless of the sterilization center's location, thereby ensuring that no streets were omitted during any campaign. Repeated coverage of the same areas ensured that dogs missed during previous campaigns, including young puppies, pregnant or lactating females, sick animals, and dogs that had escaped capture, could be sterilized during subsequent visits before reproducing and contributing to population growth.

Each route was covered during three 1.5-h shifts. The first shift commenced at 05:00 h before sunrise, when heat stress was minimal and dogs were less active and easier to capture. The second shift started at approximately 13:00 h, coinciding with the feeding activities of volunteer dog feeders, when community dogs tended to congregate at specific locations to await food. The third shift began at 16:30 h, when dogs were commonly observed following office workers and students and congregating around bus stops and shops, making them easier to capture.

The number of catching days allocated to each zone depended on necessity. Repeat visits were conducted whenever large numbers of sexually intact, free-roaming dogs of reproductive age or dogs that escaped capture were encountered. Nevertheless, both teams focused on only one zone per day, and the number of dogs captured and the duration of catching efforts remained standardized. Table 2 summarizes the number of days allocated to each zone annually.

Year | Zone | Program no. | Total | Male | Female | No. of days
2021 | 1 | 1 | 35 | 17 | 18 | 1
2021 | 1 | 2 | 139 | 28 | 111 | 4
2021 | 1 | 3 | 84 | 26 | 58 | 3
2021 | 1 | 4 | 143 | 37 | 106 | 3
2021 | 1 | 5 | 51 | 13 | 38 | 2
2021 | 1 | 6 | 77 | 21 | 56 | 2
2021 | 1 | Sub total | 529 | 142 | 387 | 15
2021 | 2 | 1 | 63 | 18 | 45 | 2
2021 | 2 | 2 | 66 | 13 | 53 | 2
2021 | 2 | 3 | 154 | 28 | 126 | 4
2021 | 2 | 4 | 92 | 23 | 69 | 2
2021 | 2 | 5 | 26 | 3 | 23 | 1
2021 | 2 | 6 | 99 | 24 | 75 | 2
2021 | 2 | Sub total | 500 | 109 | 391 | 13
2021 | 3 | 1 | 224 | 55 | 169 | 6
2021 | 3 | 2 | 160 | 34 | 126 | 4
2021 | 3 | 3 | 104 | 27 | 77 | 3
2021 | 3 | 4 | 85 | 26 | 59 | 4
2021 | 3 | 5 | 46 | 11 | 35 | 2
2021 | 3 | 6 | 146 | 35 | 111 | 6
2021 | 3 | Sub total | 765 | 188 | 577 | 25
2021 | Buffer | 1 | 24 | 11 | 13 | 1
2021 | Buffer | 2 | 0 | 0 | 0 | 0
2021 | Buffer | 3 | 0 | 0 | 0 | 0
2021 | Buffer | 4 | 0 | 0 | 0 | 0
2021 | Buffer | 5 | 0 | 0 | 0 | 0
2021 | Buffer | 6 | 0 | 0 | 0 | 0
2021 | Buffer | Sub total | 24 | 11 | 13 | 1
2022 | 1 | 1 | 99 | 27 | 72 | 2
2022 | 1 | 2 | 0 | 0 | 0 | 0
2022 | 1 | 3 | 14 | 3 | 11 | 1
2022 | 1 | 4 | 68 | 24 | 44 | 2
2022 | 1 | 5 | 56 | 18 | 38 | 2
2022 | 1 | Sub total | 237 | 72 | 165 | 7
2022 | 2 | 1 | 42 | 18 | 24 | 1
2022 | 2 | 2 | 0 | 0 | 0 | 0
2022 | 2 | 3 | 36 | 13 | 23 | 1
2022 | 2 | 4 | 89 | 28 | 61 | 3
2022 | 2 | 5 | 53 | 13 | 40 | 1
2022 | 2 | Sub total | 220 | 72 | 148 | 6
2022 | 3 | 1 | 157 | 43 | 114 | 4
2022 | 3 | 2 | 212 | 79 | 133 | 4
2022 | 3 | 3 | 53 | 16 | 37 | 2
2022 | 3 | 4 | 32 | 9 | 23 | 2
2022 | 3 | 5 | 0 | 0 | 0 | 0
2022 | 3 | Sub total | 454 | 147 | 307 | 12
2022 | Buffer | 1 | 0 | 0 | 0 | 0
2022 | Buffer | 2 | 76 | 26 | 50 | 1
2022 | Buffer | 3 | 129 | 56 | 73 | 4
2022 | Buffer | 4 | 14 | 4 | 10 | 1
2022 | Buffer | 5 | 4 | 0 | 4 | 1
2022 | Buffer | Sub total | 223 | 86 | 137 | 7
2023 | 1 | 1 | 49 | 8 | 41 | 1
2023 | 1 | 2 | 60 | 9 | 51 | 2
2023 | 1 | 3 | 52 | 15 | 37 | 2
2023 | 1 |  | 161 | 32 | 129 | 5
2023 | 2 | 1 | 71 | 19 | 52 | 2
2023 | 2 | 2 | 74 | 28 | 46 | 2
2023 | 2 | 3 | 64 | 21 | 43 | 2
2023 | 2 | Sub total | 209 | 68 | 141 | 6
2023 | 3 | 1 | 95 | 30 | 65 | 3
2023 | 3 | 2 | 62 | 11 | 47 | 2
2023 | 3 | 3 | 124 | 43 | 81 | 4
2023 | 3 | Sub total | 281 | 84 | 193 | 9
2023 | Buffer | 1 | 0 | 0 | 0 | 0
2023 | Buffer | 2 | 63 | 18 | 45 | 2
2023 | Buffer | 3 | 37 | 15 | 22 | 1
2023 | Buffer | Sub total | 100 | 33 | 67 | 3
2024 | 1 | 1 | 62 | 16 | 46 | 2
2024 | 1 | 2 | 45 | 17 | 28 | 1
2024 | 1 | 3 | 54 | 20 | 34 | 2
2024 | 1 | Sub total | 161 | 53 | 108 | 5
2024 | 2 | 1 | 78 | 18 | 60 | 2
2024 | 2 | 2 | 37 | 7 | 30 | 1
2024 | 2 | 3 | 39 | 18 | 21 | 2
2024 | 2 | Sub total | 154 | 43 | 111 | 5
2024 | 3 | 1 | 120 | 49 | 87 | 4
2024 | 3 | 2 | 81 | 12 | 69 | 2
2024 | 3 | 3 | 61 | 16 | 45 | 2
2024 | 3 | Sub total | 262 | 77 | 201 | 8
2024 | Buffer | 1 | 51 | 17 | 34 | 1
2024 | Buffer | 2 | 0 | 0 | 0 | 0
2024 | Buffer | 3 | 62 | 14 | 48 | 1
2024 | Buffer | Sub total | 113 | 31 | 82 | 2
2025 | 1 | 1 | 25 | 7 | 18 | 1
2025 | 1 | 2 | 65 | 21 | 44 | 1
2025 | 1 | 3 | 66 | 14 | 52 | 2
2025 | 1 | Sub total | 156 | 42 | 114 | 4
2025 | 2 | 1 | 50 | 7 | 43 | 1
2025 | 2 | 2 | 65 | 22 | 43 | 1
2025 | 2 | 3 | 21 | 9 | 12 | 1
2025 | 2 | Sub total | 136 | 38 | 98 | 3
2025 | 3 | 1 | 53 | 17 | 36 | 3
2025 | 3 | 2 | 73 | 23 | 50 | 2
2025 | 3 | 3 | 51 | 23 | 28 | 3
2025 | 3 | Sub total | 177 | 63 | 114 | 8
2025 | Buffer | 1 | 144 | 33 | 111 | 3
2025 | Buffer | 2 | 0 | 0 | 0 | 0
2025 | Buffer | 3 | 178 | 40 | 138 | 3
2025 | Buffer | Sub total | 322 | 73 | 249 | 6

Nylon nets mounted on stainless steel or aluminum frames were used for capturing dogs. The metallic frame had an internal diameter of 60 cm, the handle was 90 cm long, and the nylon net was 60 cm deep.

Young puppies and friendly dogs were captured manually. Catching vehicles moved slowly along predetermined routes while strictly adhering to zonal boundaries. Catchers covered narrow side streets on foot and followed the designated routes using Google Maps. Although this approach differed from conventional street surveys, its objective was to maximize capture efficiency and optimize resource utilization.

Neighboring communities were informed about upcoming clinics at least 1 week in advance through public announcements, handbills, and social media campaigns to ensure that owners who required the service did not miss the clinics due to lack of awareness or inadequate preparation. Announcements were delivered in Sinhala using loudspeakers along the intended routes, and more than 1,000 handbills were distributed before each campaign.

Data collection and record keeping

The following data were recorded for animals sterilized and vaccinated from 2020 onward: date, location, total number of dogs, total numbers of male and female dogs, numbers of owned and free-roaming dogs, numbers of male and female dogs within the owned and community dog groups, total number of cats, numbers of male and female cats, and number of ARVs administered. In addition, dogs were classified into age groups as newborn puppies with mothers, young puppies <3 months, puppies aged 3 months–1 year, sexually active dogs aged 1–6 years, and senior dogs >6 years. Dental aging techniques were used to determine age groups. Dogs found roaming on streets without an accompanying owner or visible sign of ownership, such as a collar, were initially classified as free-roaming but reclassified as owned if ownership was later confirmed.

The questionnaire survey was designed to collect information on the number of cats and dogs living in households in different regions within and outside the municipal boundary, including their age, sex, vaccination history, and sterilization status by the end of the program. Responses included information on pets brought to 49 vaccination/sterilization centers from late September to mid-October 2025, as well as other pets in the same households that were not brought to the clinics during this period. The questionnaire used for the survey is provided as Supplementary Material 1.

The survey targeted pet owners and voluntary caregivers within the Anuradhapura Municipality and surrounding Pradeshiya Sabhas covered by the project, using convenience sampling. Cochran’s formula was used to determine the minimum required sample size, with a Z-score of 1.96 for a 95% confidence interval (CI), p = 0.5 to represent maximum variability, and e = 0.05 to represent a 5% margin of error:

n = Z²p(1 − p)/e²

The total population size of the municipality and three Pradeshiya Sabhas exceeded 100,000, resulting in a minimum required sample size of 384. For the total population of 277,812 across the four divisions (Anuradhapura Municipal Council: 60,237; Manupa: 103,641; Nanupa: 70,238; and Mihinthale: 40,657), the final survey sample of 1123 was considered sufficiently representative, with an estimated 3% margin of error. Of the 1123 responses, 737 were obtained from 41 locations within the municipality and 386 from eight locations outside the municipal boundary, among 1687 pet owners who visited the clinics. The questionnaire response rate was 66.57%.

Before implementation, the questionnaire was reviewed internally for clarity and piloted among 10 team members. Responses from the pilot survey were excluded from the analysis. The authors acknowledge the sampling bias associated with surveying pet owners who visited clinics rather than conducting a door-to-door survey, because respondents were likely to represent the more responsible segment of pet owners.

Regular ARV was defined as uninterrupted annual rabies vaccination for ≥3years. Vaccination histories were verified from vaccination records provided by pet owners.

Although the survey included cat owners, data on cats were excluded from the present analysis because cat counts were limited compared with dog counts, no homeless cats were captured, and dogs were the primary focus of the project. Therefore, households with zero dogs, including households with cats only, were excluded from the analysis.

Statistical analysis

All statistical analyses were performed using RStudio version 4.4.2 (R Core Team, Vienna, Austria).

Temporal trends in dog counts: Temporal trends in counts per unit effort were quantified for nine demographic categories: total dog count; males; females; owned dogs; community dogs (free-roaming dogs cared for by the community without a specific owner); owned males; owned females; community males; and community females. Negative binomial regression was used for each category. The model count ~ Year + offset(log[number of days]) was fitted to annual totals for 2021–2025 (n = 5 years), where the offset adjusted for annual variation in the number of program days. The exponentiated Year coefficient was interpreted as the incidence rate ratio (IRR), representing the annual proportional change in catch per catching day. Each model was fitted using the glm.nb function in R v4.x. Model accuracy was assessed using the dispersion parameter θ, residual deviance/df, Pearson χ²/df, and DHARMa-scaled residuals, with consideration of uniformity, overdispersion, and zero-inflation. Influential observations were identified as Cook’s distance >0.8 (4/n = 4/5). All tests used α = 0.05.

A directed acyclic graph (DAG) was developed to justify the structure of the negative binomial regressions. The DAG identifies Year as the primary exposure and observed dog count as the outcome, with sampling effort acting as an endogenous variable driven by the underlying true population density. Including log(number of days) as an offset in each model appropriately blocks the confounding pathway associated with non-random sampling effort. Rather than adjusting for demographic categories, including sex and ownership status, as covariates in a single pooled model, analyses were stratified into nine distinct demographic categories. As illustrated in Figure 3, these demographic categories determine baseline population densities and are likely to influence catching effort. Stratification by demographic category prevents pooling bias and allows detection of effect modification, enabling calculation of distinct temporal IRRs for subsets such as community and owned dogs.

Figure 3: Directed acyclic graph showing temporal trends in dog populations stratified by demographic category.

Sex ratios: Annual male:female ratios were calculated for three populations: all dogs, owned dogs, and community dogs, to determine whether sex composition changed over time. Ratios were calculated as male dogs/female dogs, owned males/owned females, and community males/community females. Community ratios were set as missing when the denominator was zero. Ordinary least squares regression with year as the predictor was used to test temporal trends in ratios (ratio ~ year). Models were fitted separately for each category. Trends were visualized by plotting ratios against year with fitted regression lines. Analyses were conducted in R 4.4.3 using dplyr, tidyr, and ggplot2.

A DAG mapping the hypothesized underlying mechanisms (Figure 4) was developed to justify the use of unadjusted ordinary least squares regression (ratio ~ year). Improved access to free animal birth control was hypothesized to reduce the sociological tendency to abandon or reject female dogs over time. This reduction in abandonment may subsequently reduce female mortality from factors such as starvation and traffic accidents, ultimately altering the observed male:female ratio. These sociological and mortality factors act as mediators on the causal pathway between the intervention, proxied by time, and the outcome. Adjusting for these factors would introduce overadjustment bias and block the total effect. Therefore, the DAG justifies the unadjusted ordinary least squares model as the appropriate specification for capturing cumulative temporal change in sex composition.

Age-specific variation in dog counts: Age-specific variation in dogs recorded per day for each year was analyzed using negative binomial generalized linear models (GLMs) with a log link and log(number of days) as an offset, thereby accounting for variation in the number of program days each year. Separate models were fitted for four age categories: <3 months, 3 months–1 year, 1–6 years, and >6 years. Year was coded as 0 for 2021 through 4 for 2025, so that exponentiated coefficients represented annual IRRs. An IRR >1 indicated an annual increase in capture rate per day, whereas an IRR <1 indicated an annual decline. IRRs are reported with 95% CIs and p-values from Wald tests. Model fit was assessed using residual deviance divided by residual degrees of freedom, with values close to 1 indicating adequate fit. Overdispersion and zero-inflation were assessed using simulation-based diagnostics in the DHARMa package. No age-specific model showed evidence of overdispersion (all dispersion test p > 0.4) or zero-inflation (all p > 0.4). Influential observations were assessed using Cook’s distance; one to two years were influential for each model, as expected with n = 5 data points. All analyses were conducted in R v4.x using MASS, broom, and DHARMa. Significance was set at α = 0.05.

Figure 4: Directed acyclic graph showing the mechanistic pathway underlying changes in sex ratios.

A DAG mapping the operational and biological realities of the intervention (Figure 5) was constructed to justify the age-stratified negative binomial regression models. Field protocols prioritized adult dogs for sterilization over young puppies, indicating that age category directly influenced sampling effort. In addition, the biological mechanism of the sterilization program, namely prevention of new births, disproportionately affects the true population density of younger age cohorts compared with older cohorts and therefore acts as an effect modifier. Pooling all age groups into a single adjusted model would obscure these distinct dynamics because age category influences both operational effort and the biological outcome. The DAG explicitly justifies fitting separate models for each age group and using log(number of days) as an offset to account for sampling effort, thereby allowing estimation of the unique annual IRR for each age cohort.

Figure 5: Directed acyclic graph showing temporal trends by age group.

Influence of the buffer zone: To assess the influence of introducing the buffer zone on dog counts in the first three zones, dog count data were analyzed using negative binomial GLMs to account for overdispersion, with log(number of days) included as an offset. Analyses were restricted to programs that covered the zones for at least 1 day. Because programs in the early years were mostly limited to the inner zones and excluded the buffer, two periods were defined: pre-buffer intensification (2021–2022) and post-buffer intensification (2023–2025), with 2023 marking the start of intensified buffer operations. Contemporaneous buffer spillover effects in the post-2023 period were tested using the model Total ~ Zone × Period + Zone:Period:Buffer Total + offset(log[number of days]). The three-way interaction estimated the log change in zone-specific catch per additional dog sterilized from the buffer during the same program.

To test for sex-biased responses, identical models were fitted separately for male and female counts. Before/after changes were assessed using Male ~ Zone + Zone:Period + offset(log[number of days]) and Female ~ Zone + Zone:Period + offset(log[number of days]). Buffer spillover was assessed using Male ~ Zone × Period + Zone:Period:Buffer Total + offset(log[number of days]) and Female ~ Zone × Period + Zone:Period:Buffer Total + offset(log[number of days]). Buffer Total was calculated as the sum of male and female dogs sterilized from the buffer per program.

To test the hypothesis that buffer operations triggered delayed displacement into Zone 1, lagged buffer variables were created. Buffer_Run_Lag1 was coded as 1 for Zone 1 programs immediately following a buffer program, and Buffer_Gap_Lag1 represented the number of days between a buffer program and the next Zone 1 program. The models Total ~ Buffer_Run_Lag1 + offset(log[number of days]) and Total ~ Buffer_Gap_Lag1 + offset(log[number of days]) were fitted using post-2023 Zone 1 data only. The same models were fitted separately for male and female counts.

The operational dataset constrained the sample size to nine post-2023 programs per zone. Post hoc power analysis using the pwr package indicated 80% power to detect a ≥35% change in Zone 1 catch rates, given the observed dispersion parameter θ = 11.8–18.0. For buffer spillover, 80% power was available to detect IRR ≤0.992 or ≥1.008 per buffer dog, equivalent to a ≥0.8% change per dog. Thus, significant effects of IRR = 0.9943–0.9955 represented the minimum detectable effect sizes. The models were not controlled for seasonality or climate because of limited program-level metadata. Buffer Total treated all buffer removals, defined as removal of sexually active animals from the reproducing population through sterilization and not relocation, as equivalent, although spatial locations within the buffer may have varied. Model selection was based on the Akaike information criterion. Significance was assessed at α = 0.05. All analyses were performed in R 4.4.3 using MASS for negative binomial GLMs.

A comprehensive DAG (Figure 6) was constructed to justify the spatiotemporal regression models, particularly the inclusion of the Zone:Period:Buffer Total interaction, and to address limitations arising from unmeasured variables. Intensified buffer operations after 2023 were hypothesized to reduce inner-zone dog counts through two behavioral and sociological mediators: (1) reducing the availability of unwanted puppies, thereby decreasing urban dumping by pet owners, and (2) removing the mating drive among free-roaming buffer dogs, thereby preventing temporary inward migration. The DAG further illustrates that zones serve as spatial modifiers of these mediators: Zone 3 immediately borders the buffer, whereas Zone 1 is relatively insulated. Therefore, the three-way interaction captures the expected distance-decay effect of buffer spillover. Seasonality was included as an unmeasured latent node in the DAG to account for its potential role as a confounder affecting both the true population density and operational catching effort.

Figure 6: Directed acyclic graph showing the influence of the buffer zone.

Determinants of vaccination and sterilization coverage: Initially, crude vaccination and sterilization percentages were calculated for each of the 49 locations by dividing the total number of vaccinated and sterilized dogs by the total number of dogs per location. Negative binomial regression models were then fitted to counts of vaccinated and sterilized dogs, using log(total dogs) as an offset and mean dogs per household as the primary exposure to explore predictors of location-level coverage. IRRs with 95% CIs were reported. Model fit was assessed using R² from equivalent linear models of rates.

Recognizing that aggregate analysis can obscure within-location heterogeneity and is vulnerable to ecological bias, the data were refitted using binomial generalized linear mixed models (GLMMs), with household as the unit of analysis. For vaccination, the outcome was cbind(Annually Vaccinated, Not Vaccinated), and for sterilization, the outcome was cbind (Total Sterilized, Not Sterilized). Location No was included as a random intercept to account for clustering. From each model, the median odds ratio (MOR) was extracted to quantify between-location heterogeneity, and the intraclass correlation coefficient (ICC) was extracted to quantify the proportion of total variance attributable to location.

Location-specific adjusted probabilities were computed by adding each random intercept to the fixed intercept and back-transforming from the logit scale to identify priority locations. The 95% CIs were calculated for each location using the conditional variance of the random effects. Locations were ranked by adjusted sterilization probability. Priority locations were defined as those with an adjusted sterilization probability <35% and a 95% upper CI <50%. All analyses were conducted in R 4.3.2 using MASS, lme4, and merTools.

A hierarchical DAG was constructed to justify the transition from aggregate ecological analysis to GLMMs. The DAG (Figure 7) illustrates that unmeasured location-specific characteristics, referred to as latent location factors, likely act as confounders by simultaneously influencing household dog density and the probability of vaccination or sterilization. Aggregate analysis obscures these pathways and increases the risk of ecological bias. Using a binomial GLMM with Location No. as a random intercept effectively proxies these unmeasured spatial confounders. This hierarchical structure allows the model to partition variance, quantified using ICC and MOR, block confounding pathways, and calculate adjusted location-specific intervention probabilities free from aggregation bias.

Temporal patterns in TVT and skin treatments: Negative binomial GLMs, with the log of mobile clinic days as an offset, were used to test temporal trends in canine TVT and skin treatment cases reported in the municipality. This approach accounted for differences in the number of clinic days by modeling the rate of cases per clinic day.

Counts of TVT and skin treatment cases per program were used as outcome variables. Calendar year was the primary predictor and was treated as a continuous variable to estimate the annual multiplicative change in rate. Negative binomial models were selected due to overdispersion in the count data; the dispersion parameters were 4.69 for TVT and 10.55 for skin treatments.

Figure 7: Directed acyclic graph showing determinants of vaccination and sterilization coverage.

An initial attempt to include program as a random intercept using GLMM was abandoned because between-program variance was estimated to be zero, indicating no additional heterogeneity after adjustment for sampling days. Therefore, fixed-effects model results are reported. Model fit was assessed using deviance residuals and dispersion statistics.

For visualization and interpretation, annual rates were standardized to cases per 7 clinic days, calculated as sum(cases)/sum(days) × 7, to represent a typical program week and reduce stochastic variation associated with single-day rates. All analyses were performed in R v4.3.0 using the MASS package. Statistical significance was set at α = 0.05.

A DAG mapping the relationships among calendar year, sampling effort, and observed cases was formulated to justify the approach used for modeling temporal trends in TVT and skin treatments. The DAG (Figure 8) identifies Calendar Year as the primary exposure driving unobserved true disease prevalence. Because mobile clinic days may directly constrain observed counts but do not causally alter the true underlying prevalence, they were modeled as an offset, log(days), rather than as a covariate. Although unmeasured program-level factors were initially hypothesized to confound observed counts, preliminary GLMM analysis indicated zero between-program variance. Consequently, the DAG reflects a final causal structure without program-level latent confounders, justifying the use of fixed-effects negative binomial regression to estimate the annual multiplicative change in treatment rates.

Figure 8: Directed acyclic graph showing temporal trends in transmissible venereal tumor (TVT) and skin treatments.

Approximate benefit-cost analysis: Operational costs were calculated based on the operational expenditure of recent programs, including labor, transport, surgical supplies and drugs, food, and accommodation, while excluding service charges. The cost of post-exposure prophylaxis and the total cost of rabies were derived from cited sources, and the current dog population of Sri Lanka was estimated as one-eighth of the estimated human population.