Botulinum neurotoxins (BoNTs) are generated by the anaerobic bacterium Clostridium botulinum. There are seven antigenically distinct botulinum toxin serotypes (A to G), of which BoNT type A1 (BoNT/A1) has traditionally been indicated for use in cosmetics and various diseases such as axillary hyperhidrosis, chronic migraine, and neurogenic detrusor overactivity (Nigam and Nigam, 2010; Chen, 2012; Walker and Dayan, 2014; Pirazzini et al., 2017; Kwak et al., 2020).

Owing to the increased understanding of the pharmacological importance of neovascularization and neuronal signaling in maintaining the tumor microenvironment there are several studies on the anticancer effect of BoNT/A1. It has been studied in various tumor types, including glioblastoma, neuroblastoma, prostate, breast, and colorectal cancer (Nam et al., 2012; Proietti et al., 2012; Ulloa et al., 2015; Rust et al., 2016; Mittal and Jabbari, 2020). Additionally, several studies have shown that, similar to normal organs and tissues, solid tumors require innervation from the sympathetic and parasympathetic nervous systems to sustain the tumor microenvironment and facilitate cancer development, metastasis, and disease progression (Magnon et al., 2013; Coarfa et al., 2018; Faulkner et al., 2019; Reavis et al., 2020; Wang et al., 2020). Since BoNT/A1 causes denervation via cholinergic signaling pathway interference (Simpson, 2013; Kwak et al., 2020), exploring its anticancer properties is possible. In particular, various acetylcholine receptor subtypes occur in melanomas (Lucianò and Tata, 2020). More recently, it has been suggested that the presence of dorsal root ganglion neurons allows melanoma to grow significantly faster in vivo, and innervation plays a direct role in tumorigenesis by suppressing the immune response in melanoma (Keskinov et al., 2016; Reavis et al., 2020). Furthermore, a phase I clinical trial, including metastatic melanoma cancer patients, reported promising tumor inhibition activity administering a combination of pembrolizumab (PD-1 antibody) and propranolol, a non-selective beta-adrenergic blocking agent (Gandhi et al., 2021). Thus, we could conceptualize a relationship between tumor growth and innervation in melanoma. Additionally, melanoma can directly and precisely inject BoNT/A1 into the tumor site; therefore, we focused on it.

Although many studies have explored its anticancer effect and mechanism of action, quantitative exposure-response studies on BoNT/A1 have not been reported, as it is difficult to quantify the blood concentration of BoNT/A1 because of the extremely small volumes injected into the local tissue. Therefore, the regulatory agencies have approved BoNT/A1 products based on less in-depth PK data (Ravichandran et al., 2006; U.S. Food and Drug Administration, 2009; U.S. Food and Drug Administration, 2010; U.S. Food and Drug Administration, 2011). The lack of PK data is compounded by the difficulty of pharmacometrics research, including the prediction of exposure-response relationships.

To overcome these issues, the kinetic-pharmacodynamic model (K-PD), also referred to as the kinetic drug action model, can provide an alternative. The K-PD model was established based on the link between virtual pharmacokinetic (PK) and pharmacodynamic (PD) compartments. The virtual PK compartment was reversely defined, depending on the curve fit for the PD compartment observation values (Jacqmin et al., 2007; Population Approach Group Europe, 2021; Pillai et al., 2004; González-Sales et al., 2017; Lixoft, 2021a, https://mlxtran.lixoft.com/examples/k-pd-models/). The typical K-PD model structure is shown in Figure 1A. The main advantage of the K-PD model is that quantitative analysis is possible even with a lack of PK information; however, pharmacometricians can only establish the exposure-response model to determine the optimal dosing regimen.

This study used a modified Simeoni tumor growth inhibition model as a PD model for our final model build-up and we set it up based on information from the Simeoni perturbed (treated group) and unperturbed (control group) model. In the unperturbed model, tumor growth was exponential in its early phase and linear in its late phase (Simeoni et al., 2004; Mo et al., 2014; Lixoft, 2021b, https://mlxtran.lixoft.com/model-libraries/tgi-library; Li et al., 2016). Eq. 1 provides a mathematical understanding of the overall tumor growth phenomenon, including exponential and linear growth. The perturbed model combined the drug effect on tumor growth with the unperturbed model. Here, the drug effect was revealed by the drug moving from the PK compartment, i.e., the intrinsic tumor growth aspect explained by the unperturbed model could be controlled by the inhibition ability of the drug flowing into the PD compartment. Eq. 2 provides a mathematical understanding of intrinsic tumor growth and drug inhibition. d W ( t ) d t = λ 0 ⋅ W ( t ) [ 1   +   ( λ 0 λ 1 ⋅ W ( t ) ) ψ ] 1 / ψ (1) d W ( t ) d t = λ 0 ⋅ W ( t ) [ 1   +   ( λ 0 λ 1 ⋅ W ( t ) ) ψ ] 1 / ψ   − k 2   ⋅ c ( t ) ⋅ x ( t ) (2) where W(t) represents the tumor volume at time t, and λ0 and λ1 are exponential and linear growth rate constants, respectively. x(t) represents the proliferating portion of W(t). k₂ is the tumor degradation rate constant. c(t) represents the drug concentration input from the PK compartment. The power coefficient Ψ (PSI) was assumed to be 20 as a switch between the exponential and linear tumor growth phases. In Simeoni’s experience, when the value was fixed 20, it was considered to be a sharp function as a switch (Simeoni et al., 2004). The original typical structure of Simeoni’s model was shown in Figure 1B.

The first objective of this study was to develop a novel tumor growth inhibition model to explain the anticancer effect of BoNT/A1 in a melanoma syngeneic mouse model. Second, to predict the relationship of exposure-response based on our developed K-PD model, and finally, to obtain a quantitatively appropriate therapeutic dose range and regimen for single- and multiple-dose studies. Furthermore, we considered that the K-PD modeling approaches will inform pharmacometricians on the feasibility of a K-PD model build-up for BoNT/A1.

As with BoNT/A1, because the plasma drug-time-concentration profile might not always be available to correlate exposure with the biomarker or clinical endpoints during the drug discovery and development process, pharmacometricians needed to examine alternatives for cases in which PK data are unavailable. Using the syngeneic melanoma mouse model, we attempted to establish a mathematical model to explain the anticancer effect in a diseased mouse model when BoNT/A1 was intratumorally injected in a situation where PK data were not available. We attempted to combine the virtual PK part of the K-PD model with the PD part of the Simeoni TGI model. Our final K-PD model predicted the anticancer effects well. First, our final K-PD model estimated that the tumor growth rates L0 and L1 were fast for the dataset of the syngeneic melanoma mice. In addition, the virtual PK parameters were also well and stably estimated while performing curve-fitting using only the tumor volume as a dependent variable. While performing NONMEM iterations, the virtual PK parameters were also stably estimated. Meanwhile, we considered L0 and L1 omega in the final K-PD model to account for the large inter-individual variability (IIV) in tumor growth profiles (Figure 2). In developing the model, the IIV for the other parameters was considered, but it had little effect on decreasing the objective function value. Therefore, we only reflected IIV on L0 and L1 in our final K-PD model.

Based on the established K-PD model, we simulated various scenarios. First, a simulation study of a single BoNT/A1 administration was conducted. Due to the BoNT neurotoxicity in this simulation study, the dose was gradually increased to a high dose (0.0025–3200 U kg⁻¹) that could not be injected in an actual experiment, thereby securing a dose-response profile for the anti-melanoma effect—using the calculated ED₅₀ and the mouse D_max-inj. converted from the human, D_max-inj. we confirmed that the therapeutic dose range of BoNT/A1 in syngeneic melanoma mice was approximately 31.7–74 U kg⁻¹ (Figure 3; Table 2). Second, we conducted another simulation study on multiple BoNT/A1 doses from 6 to 74 U kg⁻¹ in three-dose regimens (Q.D. Q.3.D., and Q.W.). The results confirmed the dose-proportional and superior anticancer effect in the all-dose regimen scenarios (Figure 5; Table 2). Overall, we considered that the minimum meaningful dose in syngeneic melanoma mice is ≥13.8 U kg⁻¹ (0.27 U). We considered that the adequate administration interval could be Q.3.D. or Q.W. In the case of Q.D. there could be a neurotoxicological concern due to the injection of dense BoNT/A1. Meanwhile, in some BoNT/A1 products such as Dysport^® (Ipsen, United States), there was a precedent in which repeated toxicity tests were conducted by once-daily intramuscular administration for 14 days in rats, but the No Observed Adverse Effect Level (NOAEL) in this study was low such as 1 U Rat⁻¹ (∼5 U kg⁻¹) (U.S. Food and Drug Administration, 2009). Therefore, there were still concerns concerning neurotoxicity in Q.D. even though there were differences in animal species and administration sites in the tests. Furthermore, our results indicated that BoNT/A1 alone had a highly significant anticancer effect; therefore, it could be difficult to confirm the synergistic effect with other anticancer drugs such as PD-1 or PD-L1 antibody.

Fortunately, the predicted optimal dosing regimen was verified through additional in-vivo results in the B16-F10 syngeneic mouse model by performing an in vivo anticancer efficacy study with single or multiple intratumoral BoNT/A1 injections. The anticancer efficacy study conducted using BoNT/A1 injections at three doses (Supplementary Table S1) showed that most tumor volumes were within the 90% simulation interval estimated using our model (Supplementary Figure S2). Furthermore, when we compared the average observed values at 432 h, the final measurement time, with the simulation median value and interval, we found that all average observation values were within the 90% prediction interval (Supplementary Table S2). Based on the above results, we again found that our final K-PD model had appropriate predictive power. Furthermore, we applied our K-PD model to the TGI data obtained from another syngeneic model with inoculated 4T-1 breast cancer cells and validated the model, as it explained the tumor growth aspect accurately. Even when BoNT/A1’s inhibitory effect in the 4T-1 syngeneic model was slightly lower than that of B16-F10, our K-PD model exhibited a good curve-fitting for the 4T-1 tumor volume-time curve (in-house results). In addition, in a syngeneic model inoculated with MBT2 cells, which are ovarian cancer cells, the tumor growth rates, such as L1 and L0, were high and similar to that of B16-F10. Hence, for such tumor cell line types, we presume that growth rate constants similar to that of B16-F10 can be used (in-house results).

A previous study had identified the BoNT/A1 antitumoral effect in prostate cancer via the mouse xenograft experiment and clinical trials (Coarfa et al., 2018). As the result of the study, there was a reduction of tumor incidence and size where BoNT/A1’s dose was 0.45 U via intratumoral injection route to the xenograft mouse and there were meaningful findings which were not only nerve density decrease but apoptosis increases in prostate cancer when 100 U of BoNT/A1 was treated to the patients in clinical trials. Although our targeted cancer is different from this previous study, we considered that the results could be indirectly applied to our case because we have obtained a similar efficacious dose of BoNT/A1. In detail, we could expect a similar or lower effective dose (≤100 U) in the clinical study for melanoma patients considering the effective dose (100 U) of BoNT/A1 for prostate cancer patients since the minimum efficacious dose (0.27 U) in the syngeneic melanoma mouse was similar or a little lower than 0.45 U in the xenograft prostate mouse.

We are currently trying to expand the indications of BoNT/A1 by targeting melanoma that can be locally injectable. Although BoNT/A1 rarely elicits adverse events such as facial and palpebral edema, injection site pain, eye pain, erythema, psoriasis, skin infections, vertigo, nausea, fever, blepharitis, xerostomia, itching, and asthenia, we speculate that the minimum efficacious dose (0.27 U = 13.8 U Kg⁻¹ ≈ 1.1 U Kg⁻¹ for human) that we arrived at could be used for melanoma treatment in clinical trials considering the risk-to-benefit aspect.

There were some potential limitations in this study. First, our concept for the model development was applied only to B16-F10 melanoma tumors. Hence, it was unclear whether our concept for the combination model build-up could be applied to other tumors. Second, in this study, we used only tumor volume as a PD marker. However, if we had obtained the survival rate or the level of other biomarkers related to cancer growth, we would have succeeded in developing better models, such as disease progression models. Since a study to obtain these factors is planned, in this study, we focused on the development of models that could satisfactorily explain the change in tumor volume. Third, regarding the pharmacological mechanism of the anticancer effect of BoNT/A1, there have been various reports of BoNT/A1’s effects on apoptosis, phosphorylation process, and neuron denervation (Proietti et al., 2012; Bandala et al., 2013; Coarfa et al., 2018). Therefore, we have planned on performing in-vitro studies clarifying the pharmacological mechanism behind the anticancer effect. We believe that this investigation on the mechanism of tumor treatment could facilitate the development of more elaborate pharmacometrics models.

We developed a mathematical model to explain the tumor-suppressing effect of BoNT/A1 in B16-F10 melanoma cancer cells. The virtual PK compartment of the K-PD model was used as the model setup for the PK part in a situation where PK data were not available, and the PD part of the Simeoni TGI model was used to explain the tumor growth aspect in our dataset. Hence, novel tumor growth inhibition models for BoNT/A1 were established by combining the PK and PD parts of each model. Furthermore, our final K-PD models were well curve-fitted for each animal experiment data set and showed adequate predictive power for both single and repeated BoNT/A1 injection data. Moreover, we hypothesize that this approach could be scalable for pharmacometricians to use in PK-PD modeling for BoNT/A1.