A multivariate statistical analysis of the effects of styrene maleic acid encapsulated RL71 in a xenograft model of triple negative breast cancer

Multivariate statistical analyses (MVAs) have been widely applied to cancer genomics and proteomics in humans but rarely in the context of experimental cancer studies in animal models. In this study, we sought to investigate their value in the context of understanding the tumor suppressive actions of a 2^nd-generation curcumin (diferuloylmethane) analogue, 3,5-bis(3,4,5-trimethoxybenzylidene)-1-methylpiperidine-4-one (RL71) [1-9], encapsulated in styrene maleic acid (SMA) micelles (SMA-RL71; [10-13]), in animals expressing a xenograft model of triple negative breast cancer (TNBC). In this model of TNBC, SMA-RL71 (10 mg/kg, iv.) was previously shown to decrease tumor growth by 67% and modulated the expression of EGFR, Akt, mTOR, and 4EBP1 [13]. However, only univariate statistical analyses were conducted on this data set. The mechanism of action of SMA-RL71 is more likely to involve complex interactions amongst various pathways. Therefore, in this current study, the range of proteins examined was extended and multivariate statistical and data mining analyses were used, in order to identify a network of signaling proteins involved in tumor suppression [14]. We employed a combination of multiple linear regression (MLR), linear discriminant analysis (LDA), random forest classification (RFC) and cluster analyses (CA) to achieve this aim. The specific objective of using LDA and RFC was to determine whether a combination of measurements of tumor-related independent variables could be used to correctly classify animals as having received drug treatment or not. This type of classification analysis is useful in the development of biomarkers for cancer and drug-responsiveness in the treatment of cancer [15]. The specific objectives of using MLR and CA were to determine whether a combination of measurements of similar independent variables could be used to predict the value of a continuous variable, such as tumor growth, in the case MLR, or to reveal the association between different continuous variables, in the case of CA. Again, such multivariate statistical methods have been used in the field of clinical cancer, but rarely in animal studies of cancer.

The MVAs reported here are based on data partly reported previously [13]. Below is an abbreviated description of those methods. Analyses were based on n = 11 mice in the vehicle control group and n = 11 mice in the drug-treated group.

SMA-RL71 micelles were prepared as described previously [11]. SMA was used as a vehicle control by dissolving in NaOH and adjusting the pH to 7.4.

Female SCID mice (7–8 weeks old, 8/group) were inoculated s.c. into the rear flank with MDA-MB-231 cells (1 × 10⁶/0.1 ml Matrigel 50%). Once tumors reached 100 mm³, the mice were randomly allocated into treatment groups. The mice received SMA-RL71 (10 mg/kg, iv) or SMA control twice a week for 3 weeks via the tail vein. Two independent measurements of tumor volume were performed bi-weekly using electronic calipers. The mice were euthanized 24 d after treatment began and full necropsies were performed.

Tissue sections were analyzed for both microvessel density (MVD) via CD105 staining and apoptosis via the ApopTag kit as previously described [13]. Briefly, tumors were embedded in cryomatrix, sectioned (6 µm), and fixed in acetone. When slide preparation was complete, the slides were scanned with an Aperio Image ScanScope System (Leica, Chicago, IL) and analyzed by an individual who was blinded to the treatment groups. The microvessel analysis algorithm was used to quantify the MVD at a dark- and light-staining threshold of 185 and 210, respectively. The nuclear image analysis algorithm was used to quantify apoptotic stained cells as the percentage of positively stained nuclei.

To add more proteins to the data set, proliferation was quantified by determining the number of cells with positive Ki67 nuclear staining. Sections were pre-treated with antigen retrieval solution (10 mM citrate buffer with 0.05% Tween 20, pH 6.0) for 20 min at 95°C in a pre-heated jar after blocking endogenous peroxidases. Sections were then incubated with the blocking buffer in a humidified chamber for 1 h and stained with a monoclonal mouse anti-human Ki67 antibody (1:100) containing biotin, overnight at 4°C in a humidified chamber. Sections were then treated with polyclonal goat anti-mouse IgG (11 mg/L) secondary antibody for 30 min at room temperature. Negative controls were generated by substituting antigens with PBS. Sections were counterstained with haematoxylin QS, dehydrated, and DPX mounting medium was used to mount cover slips. The nuclear image analysis algorithm of the system was used to quantify the percentage (Pi) and classify the intensity (i) of positively stained proliferative nuclei and expressed as HScore = ∑Pi (i+1) [16].

Plasma thymidine kinase 1 (Tk1) levels in mice were measured by the enhanced chemiluminescent dot blot assay. Three µl of serum from control and SMA-RL7 treated mice as well as the recombinant human TK1 protein (rHTK1) standard (0.00056–0.18 µM), were applied to a nitrocellulose membrane and allowed to air dry. Membranes were then blocked with 10% non-fat milk in TBS for 1 h and washed 3× with TBS for 5 min. Membranes were incubated with anti-TK1 monoclonal antibody (1:500 in 5% BSA) at room temperature for 2 h. Membranes were then incubated (1 h) with biotinylated anti-mouse secondary antibody (1:1000) at room temperature, washed with TBST (5×) followed by streptavidin-conjugated HRP for 30 min and washed 5× with TBST. Membranes were then developed with SuperSignal substrate. Signal intensity was visualized on radiographic film and quantified with a GS-710 densitometer (Bio-Rad). PlasmaTK1 was determined using linear regression as a function of intensity and concentration of rHTK1. Protein extracts from tumors were prepared as previously described [13,16]. The density of each band was normalized to a β-tubulin loading control.

MLR: Backward MLR was used to predict tumor volume. The validity of the MLR was based on an adjusted R², which indicated the strength of the prediction [17-19]. In backward MLR, independent variables (IVs) are removed one at a time, in descending order of significance, to determine how the adjusted R² changes. Backward MLR is preferred over some other methods because it allows for the examination of the interaction between variables. Therefore, the interaction between IVs was investigated. MLR is prone to artifacts, one being that the adjusted R² increases with the number of IVs, even if they are not meaningful. The use of an ‘adjusted R^2’ partly controls for this, but nonetheless, too many IVs and excessive correlation amongst them can lead to overfitting and multicollinearity [14,20-23]. Multicollinearity was tested using a variance inflation factor and autocorrelation using the Durban-Watson statistic [14,20-23].

Because MLR is part of the general linear model (GLM), it does assume that the data conform to multivariate normality, homogeneity of the covariance matrices and are independent. Diagnostic plots in the R program were used to assess these assumptions. The overall significance of the MLR was tested using the regression ANOVA, and the significance of the individual predictors was tested using t tests [17-19,22]. All MLR analyses were performed using the programs SPSS 25 and R (version 3.6.1).

LDA: LDA is used to predict categorical variables, and was used here to predict whether the animals had received drug treatment. LDA is also part of the GLM, and therefore assumes that the data are normally distributed and are independent [14]. A form of LDA was used in which the IVs were entered together, since stepwise LDA is prone to artifacts [24]. The analysis yielded a specific LDA and a standardized canonical discriminant function that indicated which IVs are important in their relationship to the dependent variable. The statistical significance of the LDA was tested using Wilk’s λ and its validity was tested using cross-validation. Cross-validation for the LDA in this study was conducted using a leave-one-out (LOO) procedure. Simulation studies by Zavorka and Perret [25] suggest that, with k = 4 predictor variables, as was the case here, and low-moderate bivariate correlation, sample sizes in the range of n₁ = n₂ = 7−14 are sufficient (see also Lachenbruch [26]).

RFC: Decision trees have been routinely used in data mining [27-29]. Their underlying principle is that a flow-like series of questions is applied to each variable, subdividing the sample into groups that have maximal similarity by minimizing the within group variance. In the case of random forests, hundreds or thousands of decision trees can be generated simultaneously and their results combined to increase the precision of the prediction. RFC can predict a categorical variable, such as whether animals received drug treatment, from other IVs, where each tree provides a “vote” for the categorical membership, and the majority vote “wins”. RFC is largely assumption-free, and therefore it is a powerful alternative to LDA when the GLM assumptions are violated [27-30]. Furthermore, simulation studies have demonstrated that RFC can provide robust classification even with high dimensional data with small sample sizes [31]. We used RFC with 500 trees generated using 4 variables at each split. Although Breiman et al. [27] have suggested that RFC does not overfit, this view has been challenged (e.g.,[32]). We tried to avoid over-fitting by increasing the number of trees to 500 as well as by optimizing “m”, the number of variables at each split, by testing a range of values for “m”. We chose m = 4, the √p [29], as the optimal value. We used cross-validation by splitting the data into training and test data sets (70:30) and calculated out-of-bag (OOB) error based on observations that were excluded from a subset of the training data (the ‘bag’) used to produce the decision trees, and also a classification matrix in which the model based on the training data was used to predict group membership in the test data set, blind to their actual membership. The test error was not greater than the training error, which further convinced us that our model did not overfit. All RFC analyses were performed using R and the R package, Rattle [30]. In R (version 3.6.1), RFC analyses provide graphs of error as a function of the number of trees, a list of predictor variables in order or importance, as well as OOB receiver-operating characteristic (ROC) curves.

CA is a type of non-parametric analysis that is used to determine associations between variables, with no pre-determined dependent variable [14]. It does not make assumptions about the normality of the data or homogeneity of the covariance matrices; CA uses measures of ‘distance’ between variables in order to group them according to their degree of association, which is usually shown on a ‘dendrogram’ in which similarity increases as an inverse function of the y axis value [14]. The squared Euclidean distance was used with a hierarchical, Ward minimal linkage algorithm, in which clusters are formed based on the minimization of variance [14]. The data were transformed into z scores first, in order to minimize the effects of differences in scales of measurement for the different IVs.

Previous work showed that SMA-RL71 decreased MVD, as shown by CD105 staining, and also increased Apoptag-TUNEL staining and cleaved caspase-3 protein expression in tumors from treated mice [13]. However, immunohistochemical (IHC) staining of tumor sections with Ki67, a nuclear antigen expressed in proliferating cells, showed no significant difference compared to SMA controls (Fig. 1A and 1B). This was confirmed following examination of TK1 levels in the plasma, where there was no difference between the treatment groups (Fig. 1C and 1D). Thus, univariate analysis showed no difference in cell proliferation in the tumor between SMA-RL71 treatment and controls, even though tumor volume was decreased 68% [13].

While we had previously reported the role of EGFR, Akt, mTOR, and 4EBP1 proteins [13], as individually contributing to drug-mediated tumor suppression, more proteins were examined in tumors from treated mice in order to conduct MVA. The results showed that treatment with SMA-RL71 had no significant effect on the expression of Wnt between the treatment groups (Fig. 2A and 2B) but caused a 60% decrease in the expression of β-catenin (Fig. 2A and 2C). Further investigation showed that SMA-RL71 caused a 78% decrease in the expression of PKC-α in tumors (Fig. 2A and 2D), but had no significant effect on CaD1 (Fig. 2A and 2E), and PP2AA (Fig. 2A, 2F and 2G) compared to vehicle control.

The MLR adjusted R² for the prediction of tumor growth was 0.896, which was significant according to an ANOVA (F(13,8) = 14.98, P ≤ 0.0001). The Durban-Watson statistic was 2.12, which indicated a lack of autocorrelation [23]. Multicollinearity becomes a concern when the tolerance value is < 0.1 [23]. Significant predictor variables that had a tolerance ≥ 0.1, and were therefore not likely to generate multicollinearity, were drug treatment (P ≤ 0.001), Ki67 (P ≤ 0.009), PKC-α (P ≤ 0.0001), PP2AA-α (P ≤ 0.0001), PP2AA-β (P ≤ 0.0001) and CaD1 (P ≤ 0.001). The validity of this regression model was checked independently using best subsets regression and the results were similar, with R² values in the range of 86%–90% with multi-collinearity controlled for using the Mallow’s Cp index.

An LDA with a canonical correlation value of 0.88 was obtained, which was significant (Wilk’s λ (4,17) = 14.86, P ≤ 0.0001). The standardized canonical discriminant function coefficients were tumor volume 0.56, pAkt/Akt, 0.83, pEGRF, 0.71 and β-catenin, 0.57. Cross-validation demonstrated that this linear discriminant function was100% successful in classifying the animals to the correct treatment group (Table 1). The analysis was repeated using stepwise regression and similar results were obtained.

RFC was used as an alternative to LDA to predict the classification of the animals to the correct treatment groups. Five hundred trees with 4 variables at each split were used. The OOB estimate of error was only 6.25% and the error became reasonably stable after the first 400 trees. The ROC curve showed a good hit versus false alarm rate with an area under the curve (AUC) value of 0.95 (Fig. 3). In terms of variable importance, the most important variables for classifying the animals into treatment groups were pEGFR, PKC-α, PP2AA-β, pAkt, pAkt/Akt, Apoptag, β-catenin and tumor volume (Fig. 4).

Hierarchical cluster analysis showed that pEGFR clustered closely with tumor volume and Ki67, pAkt with CD105, NF-κB with PKC-α and β-catenin, and PP2AA-α with PP2AA-β, EGFR and Wnt5ab. The two largest clusters were of pEGFR, tumor volume, Ki67, pAkt and CD105 on the one hand, and Apoptag, CaD1, NF-κB, PKC-α, β-catenin, PP2AA-α, PP2AA-β, EGFR, Wnt5ab and Akt, on the other (Fig. 5).

Most animal studies investigating the role of cancer-related proteins only perform univariate analyses. As useful as this can be, differences in related proteins which function as a network, can be missed. In this study, we adopted a different approach, using MVAs to identify a network of proteins critical to growth suppression of TNBC tumors.

Using MLR, the growth of the tumor itself was strongly predicted by whether the animals had received treatment with SMA-RL71, as well as the expression of Ki67, PKC-α, PP2AA-α, PP2AA-β and CaD1. The adjusted R² of 0.896 indicated that these variables accounted for almost 90% of the variation in the growth of the tumor. The prediction potential of this group of proteins was interesting given that Ki67 (P = 0.4175), CaD1 (P = 0.2178) PP2AA-α (P = 0.1295) and PP2AA-β (P = 0.3901) were not significantly changed when analyzed by univariate methods.

It is important to note that although studies have established a relationship between Ki67 and overall survival (OS) and disease-free survival (DFS) [33-35], the predictive value of Ki67 as an indicator of chemotherapy benefit is controversial. For example, in the BR9601 adjuvant chemotherapy trial, 8 cycles of CMF (cyclophosphamide, methotrexate and fluorouracil) or epirubicin-CMF (4 cycles of epirubicin, followed by four cycles of the same CMF regimen) were administered to patients every 21 d. The primary end point was relapse-free survival (RFS) (HR = 0.36) and OS (HR = 0.30). However, Ki67 levels did not significantly interact with OS (P = 0.247) and RFS (P = 0.736) and therefore did not predict treatment outcome [36]. Similarly, in a study by the International Breast Cancer Study Group (IBCSG) IX, patients were randomly assigned to receive three 28-day courses of adjuvant CMF chemotherapy (cyclophosphamide on days 1–14, methotrexate on days 1 and 8; and 5-fluorouracil on days 1 and 8), followed by tamoxifen (57 months) or tamoxifen alone (60 months) [37,38]. CMF-tamoxifen resulted in a 5 year DFS and OS (87% and 89%) while tamoxifen also resulted in a 5 year DFS and OS (69% and 81%) [37-38]. However, there was no interaction between Ki67 levels and the response to the treatment [36,38]. Thus, Ki67 as an independent prognostic factor does not always reflect the response to treatment. However, our MVA results show that they form a biological network of molecular targets in the prediction of treatment outcome.

The role of CaD1 as a cancer metastasis associated protein was confirmed by Hou et al. [39]. Notably, CaD1 was less expressed based on IHC staining intensity and Western blotting in metastasis gastric (MKN7 and AZ521) cells compared to primary cancer cells (AGS and FU97) [39]. Additionally, IHC staining of gastric cancer tissues supported the in vitro expression of CaD1 and showed a significant decrease in CaD1 staining in lymph node metastasized tissues compared to the primary gastric tumors [39]. Furthermore, siRNA knockdown of CaD1 in AGS cells elicited an ~150% and ~50% increase in migration and invasion, respectively. In contrast, overexpression of CaD1 in AZ521 resulted in a ~50% decrease in both migration and invasion [39]. Similarly, low expression levels are found in colon cancer HCA7 cells and human breast cancer MB435S cells and a 4-fold and 7-fold increased invasion following siRNA-CaD1 treated cells compared to control cells was observed [40]. Thus, it is possible that if other studies had performed MVA, CaD1 would likely have contributed to a tumor suppressor network.

Using LDA to predict which animals had received SMA-RL71 treatment, we found a linear discriminant function that was statistically significant and 100.00% successful in correctly classifying the animals. The most important variables appeared to be tumor volume, the ratio of pAkt/Akt, pEGFR and β-catenin. A different form of classification analysis, RFC, found that pEGFR, PKC-α, PP2AA-β, pAkt, pAkt/Akt, Apoptag, β-catenin and tumor volume, were the most important predictor variables of the treatment group, although it is notable that pEGRF, pAkt/Akt and tumor volume were common to these two forms of analysis. For the RFC, the OOB estimate of error was only 6.25% and the AUC for the ROC curve, showing the hit versus false alarm rate, was 0.95. Thus, the expression of these signaling proteins was able to accurately distinguish tumor bearing mice receiving drug from those receiving vehicle with a specificity and sensitivity of 95%.

PKC-α may also be an important therapeutic target for SMARL71 in TNBC, as it was identified in the prediction of the SMA-treatment outcome during interaction with other signaling proteins in the MVAs and when analyzed as a single protein via univariate analysis. Inhibition of PKC phosphorylates Wnt/β-catenin signaling through direct phosphorylation of β-catenin at Ser45 and promotes β-catenin degradation [41]. In the present study, PKC-α was decreased by 78% following SMA-RL71 treatment in vivo. Similarly, the curcumin analogue, J1, inhibited the phosphorylation of PKC-theta in MDA-MB-231 and MCF7 cells by approximately 89% and 91%, respectively after 12 h [42]. Studies have also shown that the presence of heat shock protein 105 recruits PP2A, which in turn prevents the phosphorylation of the β-catenin degradation complex and subsequently activates Wnt-signaling that leads to cell proliferation and inhibition of apoptosis [43,44]. Shieh et al. [45] reported a time-dependent suppression of PP2A signaling by demethoxycurcumin in MDA-MB-231 cells with complete suppression at 48 h. Similarly, curcumin decreased the expression of PP2A by 50% in human rhabdomyosarcoma after 24 h, which led to an activation of mitogen-activated protein kinases and death in tumor cells [43]. Thus, there is evidence to suggest that the protein network identified does work in concert to modulate tumor growth.

CA are ‘unsupervised’ (i.e., there is no specific dependent variable to be predicted) but they explore the natural groupings of the signaling proteins. In a study by Dieninger et al. [46], 28 statistically significant m/z species were differentially expressed from 15 metastasized and 17 non-metastasized Barrett’s adenoma cases. Hierarchical clustering was performed on 10 of the most significant m/z species of the tumor tissues, which distinguished Barrett’s adenocarcinoma with lymph node metastasis from those without lymph node metastasis with 81% accuracy, following a specificity and sensitivity of 77% and 94%, respectively [42]. Also, hierarchical clustering of mass spectra of proteins identified by matrix-assisted laser desorption ionization imaging correctly separated gastric cancer and non-neoplastic mucosa [42]. CA in our studies showed that pEGFR clustered closely with tumor volume, which is consistent with the results of the LDA and RFC analyses. Also, Ki67 (anti-proliferation) clustered closely with tumor volume, which also confirmed its important role in tumor suppression. Overall, pAKT clustered with CD105 (anti-angiogenesis) while CaD1 clustered closely with Apoptag (pro-apoptosis).

We were concerned to ensure that the MVAs were valid and were not undermined by the relatively small sample size (n = 22), multicollinearity or overfitting. For the MLRs, the Durban-Watson statistic was 2.12, which indicated a lack of autocorrelation [23]. Multicollinearity becomes a concern when the tolerance value is < 0.1 [23]; however, all of the significant variables had a tolerance ≥ 0.1. We tested the validity of the MLR by repeating the analysis using a best subsets regression, which generated similar results. The validity of the LDA was tested using a ‘leave-one-out’ (LOO) cross-validation procedure, which makes over-fitting less likely, and simulation studies suggest that the sample sizes should have been sufficient given the number of predictor variables [25,26]. We tested the validity of the LDA by repeating the analysis using stepwise regression and similar results were obtained. For RFC, simulation studies have demonstrated that RFC can provide robust classification with small sample sizes [31]. We tried to avoid over-fitting by increasing the number of trees to 500 as well as by optimizing ‘m’, the number of variables at each split, by testing a range of values for ‘m’. Finally, we chose m = 4, the √p, which is recommended [29].

Statistical analysis is an integral part of drug efficacy studies, but research investigating the role of signaling pathways usually reports analysis of one protein at a time. Given the complexity of the cell signaling network in cancer, changes in one protein may not accurately reflect treatment outcome. Additionally, the interactions between signaling proteins may be missed. In this study, the MVAs have identified biological networks that drive tumor growth suppression mediated by SMA-RL71. These results should encourage others to include MVA as part of their data analysis when using in vivo cancer models.

The authors would like to thank M. Nimick and S. Taurin for technical assistance. We would also like to thank Ms. Lucy Thomsen for preparing the graphical abstract. The work was funded by a grant from the New Zealand Breast Cancer Foundation (RJR/KG) and a PhD scholarship from the University of Otago (ONKM).

The data that support the findings of this study are available from the corresponding author upon reasonable request.