We developed the protocol for this systematic review and NMA in accordance with the Preferred Reporting Items for Systematic Review and Meta-analysis (PRISMA) statement [15]. The review protocol was registered at the International Prospective Register of Systematic Reviews (registration number: CRD42021225666; accessible at https://www.crd.york.ac.uk/prospero/display_record.php?RecordID=225666) on October 1, 2020.

This systematic review and NMA of pharmacological and non-pharmacological strategies for preventing PHN was performed according to the protocol recommended by the Cochrane Collaboration [16] and reported according to the PRISMA extension for NMA guidelines [17].

We searched MEDLINE, EMBASE, Cochrane Central Register of Controlled Trials (CENTRAL), and Google Scholar from database establishment (MEDLINE; 1,946, EMBASE; 1,966) to December 2020 using search terms related to the pharmacological and non-pharmacological strategies for preventing PHN and updated it in June 2021. The search terms used for MEDLINE and EMBASE are presented in the Appendix. Two investigators (JHK and GJC) screened the titles and abstracts of the retrieved articles. The reference lists were imported to Endnote software 9.1 (Clarivate, London, UK), and duplicate articles were removed. Additional relevant articles were identified by scanning the reference lists of the articles obtained from the original search.

We only included prospective RCTs that compared two or more pharmacological and non-pharmacological strategies for preventing PHN.

The PICO-SD information is as follows:

When the data for 1, 3, and 6 months was not presented, we included the data from the nearest time point, if possible.

The following types of studies were excluded:

No language or date restrictions were applied in our study.

Two investigators (JHK and GJC) independently screened the titles and abstracts of the searched studies to identify trials that met the inclusion criteria outlined above. For articles determined to be eligible for the analysis based on their titles and/or abstracts, the full paper was retrieved. Potentially relevant studies chosen by at least one investigator were retrieved, and the full text was evaluated. Full-text articles were assessed separately by two investigators (JHK and GJC), and any disagreements were resolved through discussion. In cases where agreement could not be reached, the dispute was resolved with the help of a third investigator (HK). To minimize data duplication because of multiple reports, articles from the same author, organization, or country were compared.

The degree of agreement for study selection between the two independent investigators was computed using kappa statistics to measure the difference between the observed and expected agreements, that is, whether they were random or by chance. Kappa values were interpreted as follows: 1) less than 0: less than chance agreement; 2) 0.01-0.20: slight agreement; 3) 0.21-0.40: fair agreement; 4) 0.41-0.60: moderate agreement; 5) 0.61-0.80: substantial agreement; and 6) 0.8-0.99: almost perfect agreement [18].

Using a standardized extraction form, the following data were extracted independently by two investigators (JHK and MKK): 1) title; 2) name of the first author; 3) name of the journal; 4) year of publication; 5) study design; 6) type of pharmacological and non-pharmacological strategies; 7) dose of pharmacological agents; 8) country; 9) risk of bias; 10) inclusion criteria; 11) exclusion criteria; 12) age; 13) sex; 14) number of subjects; 15) incidence of PHN at 1, 3, and 6 months after acute HZ; and 16) pain score measured at 1, 3, and 6 months after acute HZ.

If the information were inadequate, attempts were made to contact the study authors, and additional information was requested. If unsuccessful, missing information was calculated from the available data, if possible, or was extracted from the figure using the open source software Plot Digitizer (version 2.6.8; http://plotdigitizer.sourceforge.net).

The reference lists were divided into two portions. Two investigators (HYS and MKK) completed the data extraction, one for each half of the reference list. Data extraction forms were cross-checked to verify the accuracy and consistency of the extracted data.

Risk of bias was assessed using the revised Cochrane risk of bias tool for randomized trials (RoB 2.0 version) (August 22, 2019) by two independent authors (BGK and HK) [19]. The RoB 2.0 is structured into five domains: D1) bias arising from the randomization process, D2) bias due to deviations from the intended interventions, D3) bias due to missing outcome data, D4) bias in outcome measurements, D5) bias in the selection of reported results. The overall risk of bias was evaluated. It was judged as low risk when the risk of bias in all domains was low, high when the risk of bias in at least one domain was high or the risk of bias in multiple domains was of some concern, and some concern if the overall judgement neither low nor high.

Ad-hoc tables were created to summarize the data from the included studies by showing their key characteristics and any important questions related to the review objectives. After extracting the data, the investigators determined the feasibility of the meta-analysis.

A multiple treatment comparison NMA is a meta-analysis generalization method that includes both direct and indirect comparisons of treatments. Both frequentist and Bayesian random-effects NMAs were conducted. A frequentist random-effects NMA was performed using the STATA software (version 15; StataCorp LP, College Station, TX) based on mvmeta with NMA graphical tools developed by Chaimani et al. [20].

Before conducting the NMA, we evaluated the transitivity assumption by examining the comparability of demographic data (age and sex), type of strategies for preventing PHN, and the risk of bias (low vs. removing low risk of bias from the overall risk of bias) as potential treatment-effect modifiers across comparisons.

A network plot linking the included pharmacological and non-pharmacological strategies for preventing PHN and their combination with other pharmacological and non-pharmacological strategies was constructed to indicate the types of pharmacological and non-pharmacological strategies, number of patients on different strategies, and the level of pair-wise comparisons. The nodes show the pharmacological and non-pharmacological strategies compared, and the edges show the available direct comparisons among the pharmacological and non-pharmacological strategies. The nodes and edges are weighed on the basis of the number of patients and the inverse of standard error of the effect.

We evaluated the consistency assumption within the entire network using the design-by-treatment interaction model. We also evaluated each closed loop in the network to evaluate the local inconsistencies between the direct and indirect effect estimates for the same comparison. For each loop, we estimated the inconsistency factor (IF) as the absolute difference between the direct and indirect estimates and the corresponding 95% confidence interval (CI) for each paired comparison in the loop [21]. When the IF value with 95% CI started at 0, it indicated that the direct and indirect evidence were consistent.

In the predictive interval (PrI) plot, the mean summary effects with CI were presented together with their PrIs to facilitate interpretation of the results, considering the magnitude of heterogeneity. PrIs provide an interval that is expected to encompass the estimate of a future study. A rankogram and cumulative ranking curve were drawn for each pharmacological and non-pharmacological strategy. Rankograms are the probabilities for treatments to assume a possible rank. We used the SUCRA values to present the hierarchy of pharmacological and non-pharmacological strategies for preventing PHN. SUCRA is a relative ranking measure that accounts for the uncertainty in the treatment order; it accounts for both the location and variance of all relative treatment effects. A higher SUCRA value indicates that the individual interventions have more positive results [22]. A comparison-adjusted funnel plot was used to assess the presence of small-study effects [23].

To test the robustness of the results from frequentist random NMA, we also conducted a Bayesian NMA using the fixed and random-effects models as well as Markov Chain Monte Carlo (MCMC) methods with the R statistical package gemtc [24]. Due to a lack of understanding on PHN, we used uninformative prior distributions automatically provided by the gemtc package. MCMC simulations were run using four chains with different initial values for an inferential 100,000 iterations after 50,000 burn-ins and thinning of 100. The convergence of fixed and random models derived from MCMC simulations was assessed using trace and density plots, and Gelman‐Rubin‐Brooks methods with a potential scale reduction factor (PSRF) of up to 1. And the comparing of the fit for fixed and random models was assessed using Dbar (posterior mean of the deviance), PD (adequate number of parameters), and DIC (deviance information criterion, the sum of Dbar and PD) statistics. We also calculated the SUCRA values from the Bayesian model and compared them with those in the frequentist model. We also performed a meta-regression analysis considering the potential influences of the risk of bias (low risk vs. some concern or high risk), sex, and mean age on each outcome.

The evidence grade was determined using the Grading of Recommendations, Assessment, Development, and Evaluation (GRADE) system, which involves a sequential assessment of the evidence quality, an assessment of the risk-benefit balance, and a subsequent judgment on the strength of the recommendations [25].

A total of 1,090 studies were found after searching MEDLINE, EMBASE, CENTRAL, and Google Scholar, with an additional 77 studies discovered after conducting a manual search, such as looking at references from included studies/reviews, and additional searching of related/cited articles in PubMed and Google Scholar. A total of 1,070 studies were retained after removing duplicates. After reviewing the titles and abstracts, 972 studies were excluded. In the first stage of the study selection process, the kappa value between two investigators was 0.764.

The remaining 98 studies were thoroughly examined, and 59 were eliminated for the following reasons: Meta-analysis [10,13,26-33], were irrelevant to outcomes of the study [34-55], unavailable outcomes [56-67], non-original data [68-71], grouping [72-77], and were nonrandomized studies [78-82]. As a result, 39 studies met the criteria for inclusion in this systematic review and meta-analysis (Fig. 1). In the second stage of the study selection process, the kappa value between two investigators was 0.834.

Table 1 summarizes the characteristics of the 39 studies that met the inclusion criteria. The pharmacological and non-pharmacological strategies used to prevent PHN consisted of antiviral agents (AV) [83-109]; a combination of antiviral agents and antiepileptics (AV + AE) [93,98-101,104,107-109]; electric nerve stimulation (ENS) [93,103]; a combination of antiviral agents and electric nerve stimulation (AV-ENS) [103]; a combination of antiviral agents, antiepileptics, and paravertebral block using local anesthetics and steroids (AV + AE-pv) [97,102]; a combination of antiviral agents and epidural block using local anesthetics and steroids (AV-eLS) [96]; a combination of antiviral agents and steroids (AV + SR) [88,110,111]; continuous epidural block using local anesthetics, steroids and epinepherine (EPI-LSE) [110]; steroids (SR) [92,112-115]; antiepileptics (AE) [92,93,115]; a combination of antiviral agents and intracutaneous injection using local anesthetics and steroids (AV + iLS) [105]; a combination of antiviral agents and subcutaneous injection using local anesthetics and steroids (AV + sLS) [106]; varicella zoster vaccine immunoglobulin (VZVIG) [116]; adrenocorticotropic hormone (ACTH) [112]; a combination of antivial agents and intracutaneous injection using methylene blue and local anesthetics (AV-iLM) [117]; a combination of antiviral agents and intracutaneous injection using local anesthetics (AV-iL) [105]; a combination of antiviral agents and continuous & intermittent epidural block using local anesthetics (AV-ecL) [118]; a combination of antiviral agents and intermittent epidural block using local anesthetics (AV-eiL) [118]; radiotherapy (RTx) [92]; a combination of antiviral agents, antiepileptics and stellate ganglion block using local anesthetics and steroids (AV + AE-sgLS) [101]; a combination of antiviral agents and stellate ganglion block using local anesthetics (AV-sgL) [95]; a combination of antiviral agents, antiepileptics, and intracutaneous injection using local anesthetics and steroids (AV + AE-iLS) [107]; a combination of antiviral agents, antiepileptics and cervical nerve root block using local anesthetics, steroids, and cobalamide (AV + AE-cLSC) [109]; stellate ganglion block using local anesthetics (sgL) [119]; and pulsed radiofrequency on the Gasserian ganglion (pRF) [120].

The risk of bias assessment is described in Table 2.

With regard to the outcomes of each specific datum, we generated a network plot (Fig. 2), inconsistency plot (Fig. 3), a predictive interval plot compared with the control (Fig. 4), expected mean ranking and the SUCRA values for each strategy (Fig. 5), the comparison-adjusted funnel plot (Fig. 6), and the SUCRA values from the Bayesian model compared with the SUCRA values from the frequentist model (Fig. 7). Figs. 2-7 present a summary of the results (A, B, C, D, E, and F in each figure correspond to the incidence of PHN at 3 months, incidence of PHN at 1 month, incidence of PHN at 6 months, pain score at 3 months, pain score at 1 month, and pain score at 6 months, respectively).

A total of 27 studies (3,136 patients) measured the incidence of PHN at 3 months after acute HZ. Of these, one study was separated from the loops [117]. NMA was performed excluding the abovementioned study. Thus, only 26 studies (3,072 patients) were analyzed. The network plot of all eligible comparisons for this endpoint is shown in Fig. 2A. Although all 17 management modalities (nodes) were connected to the network, two comparisons (AV and control [CTR]) were more directly compared with the other 15 nodes. No network inconsistency was observed [χ² (7) = 8.94, P = 0.257].

Ten closed loops were identified in the network after compariing the incidence of PHN at 3 months after acute HZ, but four loops (CTR/AV + AE/AE [93], CTR/ENS/AV-ENS [103], AV/SR/RTx [92], and SR/RTx/AE [92]) were formed only by multi-arm trials. Of the six closed loops, inconsistencies were observed in 1/10/12 (CTR/SR/AE) (Fig. 3A). EPI-LSE showed a lower incidence of PHN at 3 months after acute HZ than CTR in terms of 95% CI and PrI; moreover, AV + iLS and AV + AE-pv showed a lower incidence of PHN at 3 months after acute HZ than CTR, but only in terms of 95% CI (Fig. 4A, Supplementary Fig. 1A, Supplementary Table 1A). The insignificances in the 95% PrIs suggests that any future RCTs could change the significance of the efficacy of these comparisons. The rankograms and cumulative ranking curve showed that EPI-LSE followed by AV + sLS and AV + iLS had the lowest incidence of PHN at 3 months after acute HZ (Supplementary Figs. 1B, 1C). The expected mean rankings and SUCRA values of each intervention are presented in Fig. 5A. According to the SUCRA value, the incidence of PHN at 3 months after acute HZ was lower in the order of the EPI-LSE (97.0%), followed by AV + sLS (79.5%), AV + iLS (79.2%), and AV + AE-sgLS (76.7%).

The comparison-adjusted funnel plots showed that the funnel plots were symmetrical around the zero line, which suggested a less likely publication bias (Fig. 6A). The network diagnostics, using trace and density plots, showed that model convergence was valid in both fixed and random-effects models (Supplementary Figs. 1D, 1E). However, the Gelman‐Rubin‐Brooks methods with PSRF and DIC showed that the random-effects model is a slightly better fit for the data (Supplementary Figs. 1F, 1G, Supplementary Table 1). Heterogeneity was also decreased in random-effects model compared with fixed-effects model (Supplementary Table 1). Thus, we analyzed the data using a random-effects model.

The SUCRA values from the Bayesian model were similar to those from the frequentist model (Fig. 7A). When performing a meta-regression analysis using the risk of bias (low risk vs. some concern or high risk), age and sex did not improve the model fit or substantially decrease the heterogeneity, and were not considered statistically significant (Supplementary Table 3).

A total of 21 studies (2,612 patients) measured the incidence of PHN at 1 month after acute HZ. Of these, two studies were separated from the loops [117,118]. We performed an NMA excluding those studies. Thus, a total of 19 studies (2,492 patients) were analyzed. The network plot of all eligible comparisons for this endpoint is shown in Fig. 2B. Although all 15 management modalities (nodes) were connected to the network, two comparisons (AV and CTR) were more directly comparable than the other 13 nodes. No network inconsistencies were observed [χ² (4) = 7.97, P = 0.093].

Ten closed loops were identified in the network after comparing the incidence of PHN at 1 month after acute HZ. Four loops (CTR/AV + AE/AE [93], CTR/ENS/AV-ENS [103], CTR/SR/ACTH [112], and AV/ENS/AV-ENS [103]) were formed only by multi-arm trials. Of the six closed loops, inconsistencies were observed in 1/10/11 (CTR/SR/AE), and 1/2/3 (CTR/AV/AV + AE) (Fig. 3B). EPI-LSE, AV + iLS, AV + AE-pV, AV + sLS, and AV + SR showed a lower incidence of PHN at 1 month than CTR, only in terms of 95% CI not in PrI (Fig. 4B, Supplementary Fig. 2A, Supplementary Table 1B). The rankograms showed that EPI-LSE followed by AV + AE-pV had the lowest incidence of PHN at 1 month after acute HZ (Supplementary Fig. 2B). A cumulative ranking plot was drawn, and the SUCRA probabilities of the different interventions for PHN at 1 month after acute HZ were calculated (Supplementary Fig. 2C).

According to the SUCRA values, the incidence of PHN at 1 month after acute HZ was lower in the order of the EPI-LSE (98.7%), followed by AV + AE-pV (81.4%), AV + iLS (80.9%), and AV + sLS (77.0%) (Fig. 5B). The comparison-adjusted funnel plots showed that the funnel plots were symmetrical around the zero line, suggesting a lesser probability for publication bias (Fig. 6B).

The network diagnostics using trace and density plots showed that model convergence was valid in both fixed-effects and random-effects models (Supplementary Figs. 2D, 2E). However, the Gelman‐Rubin‐Brooks method with the PSRF and DIC showed that the random-effects model is a slightly better fit for the data (Supplementary Figs. 2F, 2G, Supplementary Table 2). Thus, we analyzed the data using a random-effects model. The SUCRA values from the Bayesian model were similar to those from the frequentist model, demonstrating the robustness of our analysis (Fig. 7B). When performing a meta-regression analysis using the risk of bias (low risk vs. some concern or high risk), age and sex did not improve the model fit or substantially decrease the heterogeneity, and were not considered statistically significant (Supplementary Table 3).

A total of 17 studies (2,502 patients) measured the incidence of PHN at 6 months after acute HZ. Although all 17 management modalities (nodes) were connected to the network, three comparisons (AV, AV + AE and CTR) were more directly comparable than the other 14 nodes (Fig. 2C). There was no evidence of network inconsistency [χ² (6) = 8.72, P = 0.190]. Thirteen closed loops were identified in the network after comparing the incidence of PHN at 6 months after acute HZ, but four loops (CTR/AV + AE/AE [93], CTR/ENS/AV-ENS [103], AV/SR/RTx [92], and SR/RTx/AE [92]) were formed only by multi-arm trials. Of the six closed loops, inconsistencies were observed in the 1/2/3 loop (CTR/AV/AV + AE) (Fig. 3C). AV + AE-pv showed a lower incidence of PHN at 6 months after acute HZ than CTR, only in terms of 95% CI (Fig. 4C, Supplementary Fig. 3A, Supplementary Table 1C). The AV + AE-cLSC showed marginal significance.

The rankograms and cumulative ranking plot showed that AV + AE-sgLS, AV + AE-cLSC, AV + AE-pv, and AV + AE-iLS had the lowest incidence of PHN at 6 months after acute HZ (Supplementary Figs. 3B, 3C).

According to the SUCRA value, the incidence of PHN at 6 months after acute HZ was lower in the order of the AV + AE-sgLS (87.8%), followed by AV + AE-cLSC (81.7%), AV + AE-pv (77.9%), and AV + AE-iLS (76.0%) (Fig. 5C). The comparison-adjusted funnel plots showed that the funnel plots were symmetrical around the zero line, suggesting a lesser probability for publication bias (Fig. 6C). The network diagnostics, using trace and density plots, showed that model convergence was valid in both fixed-effects and random-effects models (Supplementary Figs. 3D, 3E). However, the Gelman‐Rubin‐Brooks methods with PSRF and DIC showed that the random-effects model is a slightly better fit for the data (Supplementary Figs. 3F, 3G, Supplementary Table 2). Thus, we analyzed the data using a random-effects model.

The SUCRA values from the Bayesian model were similar to those from the frequentist model, demonstrating the robustness of our analysis (Fig. 7C). When performing a meta-regression analysis using the risk of bias (low risk vs. some concern or high risk), age and sex did not improve the model fit or substantially decrease the heterogeneity, and were not considered statistically significant (Supplementary Table 3).

A total of 14 studies (2,047 patients) measured the pain score at 3 months after acute HZ. Of these, three studies [104,114,119] did not report any information on the degree of scattering. Thus, we attempted to contact the study authors but could not obtain information on the degree of scattering. As two studies were separated from the loops [110,120], we performed an NMA, excluding those studies. Thus, only nine studies (1,225 patients) were analyzed. The network plot of all eligible comparisons for this endpoint is shown in (Fig. 2D). Eight management modalities (nodes) were connected to the network. One closed loop (the 1/2/5 loop) (AV/AV + AE/AV + AE-pv) was identified in the network after comparing the pain score at 3 months after acute HZ, which shows no evidence of inconsistency (Fig. 3D).

A PrI plot and league table comparing all groups were drawn (Supplementary Fig. 4A, Supplementary Table 1D). The rankograms and cumulative ranking plot showed that AV + iLS and ENS had the lowest pain score at 3 months (Supplementary Figs. 4B, 4C). According to the SUCRA values, the pain score at 3 months was lower in the order of the ENS (93.0%) followed by AV + iLS (91.0%), followed by AV + AE-sgLS (65.9%) (Fig. 4E). The comparison-adjusted funnel plots showed that the funnel plots were symmetrical around the zero line, which suggested a lesser probability for publication bias (Fig. 6D). The network diagnostics using trace and density plots showed that model convergence was valid in both fixed-effects and random-effects models (Supplementary Figs. 4D, 4E). However, the Gelman‐Rubin‐Brooks methods with PSRF and DIC showed that the random-effects model is a slightly better fit for the data (Supplementary Figs. 4F, 4G, Supplementary Table 2). Thus, we analyzed data using random-effects model.

The SUCRA values from the Bayesian model were similar to those from the frequentist model, demonstrating the robustness of our analysis (Fig. 7D). When performing a meta-regression analysis using the risk of bias (low risk vs. some concern or high risk), age and sex did not improve the model fit or substantially decrease the heterogeneity, and were not considered statistically significant (Supplementary Table 3).

A total of 21 studies (2,758 patients) measured the pain score at 1 month after acute HZ. Of these, three studies [104,114,119] did not report any information on the degree of scattering. Thus, we attempted to contact the study authors, but could not obtain information on the degree of scattering. As three studies were separated from the loops [116,117,120], we performed NMA excluding those studies. Thus, only 15 studies (2,385 patients) were analyzed. The network plot of all eligible comparisons for this endpoint is shown in (Fig. 2E). Although all 12 management modalities (nodes) were connected to the network, two comparisons (AV and AV + RE) were more directly comparable than the other 10 nodes.

There was no evidence of network inconsistency [χ² (1) = 0.14, P = 0.713]. One closed loop (the 1/2/5 loop) (AV/AV + AE/AV + AE-pv) was identified in the network after comparing the pain scores at 1 month after acute HZ, which showed inconsistencies in this loop (Fig. 3E).

A PrI plot and league table, comparing all groups, were drawn (Supplementary Figs. 5A, Supplemeantary Table 1E).

The rankograms and cumulative ranking plot showed that AV + iLS and AV + AE-sqLS had the lowest pain score at 1 month (Supplementary Figs. 5B, 5C). According to the SUCRA values, the pain score at 1 month was lower in the order of the AV + AE-sgLS (71.8%), followed by AV + iLSL (66.1%), AV-sgL (63.0%), and AV + AE-iLS (61.0%) (Fig. 5E). The comparison-adjusted funnel plots showed that the funnel plots were symmetrical around the zero line, suggesting a lesser probability for publication bias (Fig. 6E).

The network diagnostics, using trace and density plots, showed that model convergence was valid in both fixed-effects and random-effects model (Supplementary Figs. 5D, 5E). However, the Gelman‐Rubin‐Brooks methods with PSRF and DIC showed that the random-effects model is a slightly better fit for the data (Supplementary Figs. 5F, 5G, Supplementary Table 2). Thus, we analyzed data using the random-effects model.

The SUCRA values from the Bayesian model were similar to those from the frequentist model, demonstrating the robustness of our analysis (Fig. 7E).

When performing a meta-regression analysis using the risk of bias (low risk vs. some concern or high risk), age and sex did not improve the model fit or substantially decrease the heterogeneity, and were not considered statistically significant (Supplementary Table 3).

A total of 10 studies (1,407 patients) measured the pain score at 6 months after acute HZ. Of these, two studies [97,114] did not report any information on the degree of scattering. Thus, we attempted to contact the study authors but could not obtain information for the degree of scattering. As two studies were separated from the loops [110,120], we performed an NMA excluding those studies. However, as the ENS/AV/AV + iLS loop and AV + AE-pV/AV + AE/AV + AE-iLS/AV + AE-cLSC/AE-sgLS were separated, an NMA was not performed (Fig. 2F).

Five outcomes were evaluated using the GRADE system. The evidence quality was moderate for the incidence of postherpetic neuralgia at 1, 3, and 6 months after acute HZ and the pain score measured at 1 month after acute HZ, whereas the evidence quality was low for pain score measured at 3 months after acute HZ (Table 3).