We developed an in silico approach to model B16F10 melanoma cell response to a helium atmospheric pressure plasma jet (APPJ) or/and doxorubicin drug (DOX). The in silico model is informed by relevant data from previously published in vitro experiments (cancer cell viability), providing detailed information on (i) cell population number (Ncell) development during incubation and (ii) probability values for apoptosis (%PApoptosis) and mitosis (%PMitosis) following cell subjection to the plasma-conditioned RPMI-1640 medium (PC-RPMI), DOX, and DOX combined with APPJ. When treating cancer cells with PC-RPMI and DOX separately, at the smallest plasma duration (dPlasma = 15 s) and DOX concentration (cDOX = 0.05 μM), only a small decline in Ncell, an increase in %PApoptosis, or/and a decrease in %PMitosis are measured with respect to the control conditions (non-treated cancer cells). However, cell cytotoxicity is increasingly enhanced with increasing dPlasma and cDOX up to 120 s and 0.5 μM, respectively. At those highest values studied in in silico, simulated %PApoptosis are significantly larger than %PMitosis, resulting in a severe decrease in Ncell compared to control in agreement with the corresponding in vitro experiments. Furthermore, cell treatments combining the smallest two cDOX (0.05 and 0.1 μM) with dPlasma = 15 s result in smaller Ncell, larger %PApoptosis, and lower %PMitosis compared to PC-RPMI and DOX effects alone. The present in silico model is particularly useful in the plasma (cancer) medicine field since it can effectively simulate and quantify responses of various cancers to APPJ or/and cancer drugs being strongly complementary to in vitro experiments.
Atmospheric pressure plasma jets (APPJs) are non-equilibrium discharges with very promising applications in biomedicine. APPJs generate localized and controlled fluxes of reactive oxygen and nitrogen species (RONS), (V)UV radiation, and transient electric fields that can affect the functionality of biological matter.1,2 At the same time, their gas temperature remains close to room temperature,3,4 thus avoiding thermal damage of living tissues. Specifically, their efficacy and selectivity against cancer have been demonstrated in in vitro and in vivo (or in ovo) studies through: (i) plasma radiation applied on cancer cells maintained in liquid media or solid tumors in mice models (direct treatment) and (ii) buffered solutions (or cell-growth media) irradiated by APPJs with the conditioned liquid (containing newly formed RONS) being put in interaction with cancer cells or injected into solid tumors (indirect treatment).5–7 Although APPJ treatments have shown effective against different cancers, e.g., normal transformed lung fibroblasts (MRC5Vi)8 and squamous pharynx cell carcinomas (Fadu),9 they were not as effective in others [human colon cancer cells (HCT116),8 squamous tongue cell carcinomas (CAL27),9 respectively]. Thus, increasing the duration of the plasma treatment or using APPJ as adjuvant therapy with other cancer treatment strategies (chemotherapy, pulsed electric fields, etc.) may offer improved therapeutic outcomes.5,10–12 It should be noted that APPJ treatments differ significantly from traditional invasive therapies such as chemotherapy and radiotherapy. This is due to the fact that APPJ generates key RONS inducing cytotoxic effects in cancer cells without affecting the healthy cells.6
To date, most research works about APPJ effects on neoplasia consider in vitro and in vivo investigations with much fewer reporting in silico (computer models) studies. In silico works in plasma medicine have focused on the description of physical, chemical, and biological processes that determine plasma effects on cells, healthy or wounded skin, liquid-covered tissues, solid tumors, etc.13 Zero-dimensional models have been used to investigate chemical kinetics occurring, for instance, in the plasma–gas and/or plasma–liquid (inter)phases.14,15 Computational fluid dynamics models have been employed to describe plasma interactions with conductive and non-conductive surfaces and provide insights on discharge dynamics, key RONS production (such as OH, NO, NO2, H2O2, and 1Ο2), electron and electric field properties, and their transport efficiency to biological matter.16–18 Molecular dynamics simulations have been introduced to study interactions of plasma-generated RONS with biomolecules (in or around cancer cells) and DNA, RONS penetration through cell membranes, plasma oxidation of proteins, and so on.19,20 Despite the progress made thus far in in silico, no models of cancer cell dynamics induced by APPJs have yet been proposed, the only exception being the work by Murphy et al.,21 using a 3D hybrid discrete–continuum model to investigate tumor growth and apoptosis upon exposure to plasma-produced RONS. Such models can be very complementary to in vitro and in vivo experiments for quantifying and better understanding effects of APPJs (or other treatments) on various cancers in a noninvasive manner (e.g., less animal studies).
This contribution proposes an original in silico approach to quantify melanoma cancer cell response to helium APPJ or/and doxorubicin drug (DOX). The in silico framework is informed by relevant data from in vitro experiments on C57BL/6 murine B16F10 melanoma cancer cells, following the results published by Pefani-Antimisiari et al.10 Thus, APPJ and DOX features, cell treatment protocols, and viability/cytotoxicity assessments are described there. Details on the electrical and optical features of the APPJ are also given in the supplementary material. The selected in vitro experiments from Ref. 10 on which the construction of the in silico model was based are graphically represented in Fig. 1. In brief, Fig. 1(a) depicts DOX or/and APPJ applications to the RPMI-1640 medium (hereafter denoted RPMI) being in contact with the cancer cells. Following their exposure to DOX or/and APPJ, cancer cells are incubated for 48 h (t48h), and their viability is determined afterwards [Fig. 1(b)]. Based on this information, three cell-treatment scenarios, experimentally investigated in Ref. 10 (see Fig. 2 there), are also considered in the in silico model: (i) B16F10 cells in RPMI are exposed only to DOX (concentration: cDOX = 0.05–0.5 μM); (ii) RPMI without cancer cells is subjected to APPJs only (duration: dPlasma = 15–120 s) and the resulting plasma-conditioned RPMI (PC-RPMI) is put in contact with B16F10 cells (indirect plasma treatment); (iii) B16F10 cells in RPMI are exposed to DOX (cDOX = 0.05 and 0.1 μM) and then subjected to APPJ for 15 s (direct plasma treatment). In each “treatment case” scenario, the evolution of B16F10 cell number during incubation is simulated with in silico results fitted to the experimental data [like for the control in Fig. 1(c)].
(a) Graphical representation of the in vitro experimental setup of B16F10 cells treated with DOX or/and APPJ; cancer cells were maintained in a 0.5 ml of the RPMI-1640 medium. (b) Quadruplicate cell treatments were performed using the four corner-wells of different 24-well plates (S1–S4). After each treatment, cells were incubated for 48 h (t48h), and their viability was measured at t48h (MTT colorimetric assay). (c) Experimental (red), counted after t0, t4h, t24h, and t48h (Trypan blue exclusion assay, calculated from Fig. S2 in Ref.10), and simulated (black; this work) B16F10 cell number during incubation under control conditions (i.e., non-treated cells). (d) Cancer-cell mechanisms incorporated in the in silico cancer model. All in vitro experiments were performed by Pefani-Antimisiari et al., Sci. Rep. 11, 14788 (2021). Copyright 2021 The Author(s), licensed under a Creative Commons Attribution (CC BY) license.10
(a) Graphical representation of the in vitro experimental setup of B16F10 cells treated with DOX or/and APPJ; cancer cells were maintained in a 0.5 ml of the RPMI-1640 medium. (b) Quadruplicate cell treatments were performed using the four corner-wells of different 24-well plates (S1–S4). After each treatment, cells were incubated for 48 h (t48h), and their viability was measured at t48h (MTT colorimetric assay). (c) Experimental (red), counted after t0, t4h, t24h, and t48h (Trypan blue exclusion assay, calculated from Fig. S2 in Ref.10), and simulated (black; this work) B16F10 cell number during incubation under control conditions (i.e., non-treated cells). (d) Cancer-cell mechanisms incorporated in the in silico cancer model. All in vitro experiments were performed by Pefani-Antimisiari et al., Sci. Rep. 11, 14788 (2021). Copyright 2021 The Author(s), licensed under a Creative Commons Attribution (CC BY) license.10
Simulated B16F10 cell number after (a) 4 h (t4h) and (b) 48 h (t48h) of incubation. (c) Cell apoptosis (black) and mitosis (red) probability values evaluated in silico after five statistically independent ABM simulations. Black and red dashed lines represent mean apoptosis and mitosis probability values from those five simulations, respectively. This image refers to the control experiment (non-treated cancer cells).
Simulated B16F10 cell number after (a) 4 h (t4h) and (b) 48 h (t48h) of incubation. (c) Cell apoptosis (black) and mitosis (red) probability values evaluated in silico after five statistically independent ABM simulations. Black and red dashed lines represent mean apoptosis and mitosis probability values from those five simulations, respectively. This image refers to the control experiment (non-treated cancer cells).
The proposed in silico approach employs the agent-based modeling method22 (ABM) to simulate the in vitro experiments of B16F10 cancer cells performed by Pefani-Antimisiari et al.10—a relevant in silico in the concept approach to simulate glioblastoma cells is presented by de Montigny et al.23 In brief, the ABM is a complex-system method that models autonomous and interactive agents, denoted as particles positioned in space.24 As with the in vitro, a two-dimensional space is assumed where agents can freely move in a square plane (4.84 mm2), providing restrictions for agent–agent overlap and collapse. Periodic boundary conditions are applied at the sides of the 2D square, while temporal resolution in our simulations assumes a 1-h time increment [as shown in Fig. 1(c); black]. The behavior of each agent is determined after accounting in simple rules and interactions with other agents, as well as with respect to external stimuli, here PC-RPMI, DOX, and DOX combined with APPJ. Thus, agents represent here individual B16F10 cancer cells. They are modeled to grow, migrate, divide, and die [apoptosis; see Fig. 1(d)] with the latter two mechanisms of cell behavior being modulated by the effect of the cytotoxic drug, DOX, while PC-RPMI and APPJ modulate apoptosis only. Contrary to this, cell growth and migration are unaffected by those external stimuli. The likelihood (probability) for cell migration and cell growth is fixed to 99% and 50%, respectively, whereas the probability for cell apoptosis and division (either programmable or induced by DOX and/or PC-RPMI/APPJ) is parametrically investigated in the results that follow. Cell migration combines a random walk movement (average rate: 1 μm/min) and O2-gradient chemotaxis (maximum rate: 5 μm/min). The cells' mechanical interactions (when in proximity) are modeled through Newton repellant forces, while cell adhesion is neglected to reflect the conditions of matrix cells' suspension in the RPMI. The ABM describes drug and plasma-generated RONS as a “substance” whose concentration is assumed homogeneous everywhere in the simulated domain of the well; hence, their transport inside the well is almost instantaneous (diffusion is neglected) while no dissipation is accounted in the extracellular space. Regarding the in silico initial conditions, Ncell = 700 cells are considered to match the experimentally measured density of ∼25 × 103 cells per well (considering an average 25 ± 5 μm cell diameter),10 while cDOX and dPlasma are set according to reflect each “treatment case” scenario. The ABM is implemented in the C++ project invitro_neuro,25 which employs the open-source software platform BioDynaMo.26 In the supplementary material, we provide more details on the assumptions and details of the model.
Figures 2(a) and 2(b) depict indicative simulated cell population number (Ncell) developments at t4h and t48h under control conditions [a full animation (t0—t48h) is available in the supplementary material; cell colors represent different cell ages, as explained there]. Ncell is obtained from Fig. 1(c) by appropriately fitting in silico (black; this work) to in vitro data (red; calculated based on Fig. S2 of Ref. 10). Indeed, at t4h in Figs. 1(c) and 2(a), the simulated Ncell is very close to the initial cell number, as for the in vitro experiment. Then, at t24h the model predicts a ∼7-fold increase in Ncell [∼36% lower than the corresponding in vitro value; see Fig. 1(c)], while at t48h, it predicts a ∼30-fold increase in Ncell with respect to the initial cell number [Figs. 1(c) and 2(b)], being only ∼12% higher than for the experiment. Therefore, the in silico model is very complementary to the in vitro experiments. In addition, the model provides real-time quantification [1-h time increment; black data in Fig. 1(c)] and visualization [Figs. 2(a) and 2(b)] of Ncell during incubation. This is practically not feasible with in vitro experiments [red; Fig. 1(c)] due to the need of using time-consuming cell-viability assays (e.g., MTT and Trypan blue) and expensive B16F10 melanoma cell vials. Furthermore, to appropriately simulate Ncell development between t0 and t48h, the model easily encompasses and analyzes different probability values for cancer cell apoptosis (%PApoptosis) and mitosis (%PMitosis) (or other cell mechanisms if needed). Then, for each “treatment case” scenario, the parameter values of %PApoptosis and %PMitosis for which the in silico results best fit to the in vitro data [e.g., Fig. 1(c)] are obtained through trial and error. This further reduces the need of sophisticated flow cytometry experiments for exposing different cancer cell mechanisms (e.g., annexin V and propidium iodide staining for detecting apoptotic and necrotic cells, respectively).10 The values of %PApoptosis and %PMitosis for five statistically independent in silico experiments of the control melanoma cell development are presented in Fig. 2(c). %PApoptosis varies between 15.31 ± 0.14% and 18.92 ± 0.24% giving an average %PApoptosis of 17.54 ± 1.72%. However, %PMitosis obtains clearly larger values between 22.74 ± 0.08% and 24.98 ± 0.04%, resulting in an average %PMitosis of 23.72 ± 1.1%, i.e., ∼6% higher than for the apoptosis. Using those average probabilities, the calculated apoptosis–to–mitosis ratio (RA–M) is ∼0.73 ± 0.04 under control conditions. Thus, when cancer cells remain untreated, they are more prone to undergo mitotic than apoptotic mechanisms, which leads to their substantial proliferation at t48h [Fig. 2(b)]. In the following, this RA–M value will be referred to as a critical threshold to identify and assess the impact on B16F10 cell number/viability of PC-RPMI, DOX, and DOX combined with APPJ. Specifically, larger RA–M than the previous threshold would mean enhanced %PApoptosis and/or reduced %PMitosis for cancer cells and, thus, smaller viability and Ncell. In the following figures, mean values and standard deviation bars correspond to five statistically independent computational experiments of the in silico model.
Figure 3 illustrates the simulation results about the effect on B16F10 cell dynamics and mechanisms of PC–RPMI under different dPlasma. All data refer to t48h, i.e., end of incubation. Comparison of Fig. 3(a) against Fig. 2(b) (control) shows that, for dPlasma = 15 s, PC-RPMI has a negligible role in the decrease in the cell viability/number. However, an eightfold increase in dPlasma [Fig. 3(b)] drives into an approximately 20 times smaller Ncell with respect to the control experiment [also see Figs. S1(a) and S1(b) in the supplementary material]. Therefore, B16F10 cell viability is significantly reduced due to their interaction with larger amounts of long-lived RONS in PC-RPMI (such as H2O2, NO2−, NO3−). Those species have demonstrated synergism in inducing strong cytotoxic effects on cancers depending on their concentration.2,6–9 Their concentration in PC-RPMI was not measured in the present work due to the unavailability of relevant colorimetric assays. Furthermore, Fig. 3(c) illustrates the dependence of RA–M on dPlasma. At dPlasma = 15 s, RA–M is almost the same with the control threshold (dashed red), resulting in similar cell dynamics [Figs. 3(a) and 2(b)]. As so, dPlasma = 15 s cannot lead to significant modification in the medium's chemistry, and it is not suitable for enhanced B16F10 cell apoptosis and viability decrease. The PC-RPMI impact starts to become noticeable after dPlasma = 30 s. Indeed, the corresponding RA–M is 0.93 ± 0.03, and it goes up to 1.26 ± 0.03 at dPlasma = 120 s, which is approximately 27% and 73% larger, respectively, as compared to the control threshold. Therefore, by fitting the in silico to the in vitro data, simulations underpin that %PApoptosis elevates when increasing the plasma duration with a rate of 0.11% approximately [see Fig. S1(c)]. Note that in Fig. 3, %PMitosis was kept the same as that of the control [23.72% ± 1.1%; Fig. S1(c)]. In fact, we are not aware of any published study reporting on enhanced or suppressed effects of PC-RPMI on B16F10 murine melanoma cell mitosis with increasing dPlasma.
Simulated B16F10 cell number at t48h following subjection to PC-RPMI at (a) dPlasma = 15 s and (b) dPlasma = 120 s. (c) Cancer-cell apoptosis-to-mitosis ratio evaluated in silico as a function of dPlasma.
Simulated B16F10 cell number at t48h following subjection to PC-RPMI at (a) dPlasma = 15 s and (b) dPlasma = 120 s. (c) Cancer-cell apoptosis-to-mitosis ratio evaluated in silico as a function of dPlasma.
The in silico results demonstrating the DOX effect on B16F10 cell dynamics and mechanisms are illustrated in Fig. 4. The simulated Ncell corresponding to the smallest cDOX [0.05 μM; Fig. 4(a)] is ∼1.4 times smaller than that of the control [Fig. 2(b); also see Fig. S2(a) in the supplementary material]. When compared to the smallest dPlasma in Fig. 3(a), this cDOX is more effective in reducing the cell number/viability. In fact, at dPlasma = 15 s, it is expected that the amount of long-lived RONS (H2O2, NO2−, NO3−, etc.) generated by the APPJ in the RPMI is not sufficient to induce significant apoptosis in cancer cells.8,9,12 Furthermore, at the largest cDOX of 0.5 μM in Fig. 4(b), a drastic drop in Ncell is revealed [about 22-fold reduction with respect to control; Fig. S2(b)]. These results can be further analyzed by studying the dependence of RA–M on cDOX [Fig. 4(c)]. First, the RA–M corresponding to cDOX = 0.05 μM is by ∼7% larger than for the control threshold (dashed red). From Fig. S2(c), this increase is due to slightly larger %PApoptosis and smaller %PMitosis at this drug concentration. Furthermore, when cDOX is increased (0.1–0.2 μM), RA–M becomes ∼13.5% larger than for the control threshold. Again, this augmentation is justified by persistently larger and smaller %PApoptosis and %PMitosis, respectively, with increasing cDOX [Fig. S2(c)]. Finally, at cDOX = 0.5 μM, the model demonstrates a dramatic increase in RA–M, which is ∼59% larger than for the control threshold. From Fig. S2(c), cDOX = 0.5 μM exhibits the smallest %PMitosis and largest %PApoptosis among all cDOX studied, while it is the only drug concentration for which %PApoptosis is clearly larger than %PMitosis, leading to enhanced cell cytotoxicity and viability decrease. It should be noted that an RA–M value of 1.16 ± 0.05 in Fig. 4(c) (cDOX = 0.5 μM) could also be achieved when using PC-RPMI at dPlasma = 90–100 s, since for dPlasma = 60 and 120 s, the corresponding RA–M are 1.01 ± 0.02 and 1.26 ± 0.03 [Fig. 3(c)], suggesting enhanced production of long-lived RONS in PC-RPMI.
Simulated B16F10 cell number at t48h following exposure to DOX in the RPMI at (a) cDOX = 0.05 μM and (b) cDOX = 0.5 μM. (c) Cancer-cell apoptosis-to-mitosis ratio evaluated in silico as a function of cDOX in the RPMI.
Simulated B16F10 cell number at t48h following exposure to DOX in the RPMI at (a) cDOX = 0.05 μM and (b) cDOX = 0.5 μM. (c) Cancer-cell apoptosis-to-mitosis ratio evaluated in silico as a function of cDOX in the RPMI.
Figure 5 presents in silico results, which demonstrate the effect on B16F10 cell dynamics of combined DOX and direct APPJ treatment. To perform the investigation, the smallest two cDOX from Fig. 4(c) (0.05 and 0.1 μM) and dPlasma from Fig. 3(c) (15 s) were used. This was done to obtain measurable cell numbers, since combining higher cDOX and dPlasma would result in a critical decrease in cell viability.10 Therefore, it would be more difficult to discriminate the impact on cancer cells of the combined action of DOX and APPJ. Figure 5(a) illustrates simulated cell dynamics corresponding to cDOX = 0.1 μM. The simulated Ncell in this scenario is smaller than that in Fig. 4(a) (cDOX = 0.05 μM), as imposed by the higher corresponding RA–M in Fig. 4(c). Specifically, the simulated Ncell at CDOX=0.1 μM is ∼2 times smaller than that for the control [Fig. 2(b)]. However, by additionally subjecting cancer cells to dPlasma = 15 s (in the presence of DOX), see Fig. 5(b), Ncell drops by about 4.5 times with respect to the control [see Fig. S3(a) in the supplementary material). Thus, our simulations demonstrate a significant enhancement in cell cytotoxicity of 2.25 and ∼4 times, compared to that induced only by DOX or PC-RPMI, respectively. In fact, when DOX-treated cancer cells are further subjected to direct APPJ treatment, it is possible that they interact with plasma-generated short-lived RONS (O, OH, NO, 1O2, etc.), UV radiation and/or pulsed electric fields, which may also induce cytotoxic effects.1,2,9,11 Finally, the dependence of RA–M on the combined action of DOX and APPJ is depicted in Fig. 5(c). For cDOX = 0.05 μM and dPlasma = 15 s, RA–M is 0.87 ± 0.03, while it goes up to 0.91 ± 0.05 for cDOX = 0.1 μM and dPlasma = 15 s. This corresponds to an increase in ∼19% and 24.5% with respect to the control threshold (dashed red). Therefore, combining DOX with direct APPJ treatments results in a larger B16F10 cell cytotoxicity. As shown in Fig. S3(b), our model captures the effect of combined DOX with direct plasma treatment in enhancing %PApoptosis and suppressing %PMitosis of B16F10 melanoma cells. This finding agrees with recent in vitro experiments performed by Pefani-Antimisiari et al.,10 where DOX and direct APPJ treatment demonstrated a synergistic action on melanoma cell viability decline.
Simulated B16F10 cell number at t48 (a) in the presence of 0.1 μM DOX in RPMI, and (b) 0.1 μM DOX in RPMI treated for 15 s with APPJs. (c) Cancer-cell apoptosis-to-mitosis ratio evaluated in silico for two scenarios of combined treatments using DOX and APPJ.
Simulated B16F10 cell number at t48 (a) in the presence of 0.1 μM DOX in RPMI, and (b) 0.1 μM DOX in RPMI treated for 15 s with APPJs. (c) Cancer-cell apoptosis-to-mitosis ratio evaluated in silico for two scenarios of combined treatments using DOX and APPJ.
To summarize, in this work, we developed an in silico agent-based model (ABM) to simulate B16F10 murine melanoma cancer cell response to helium APPJ or/and DOX. The in silico model was informed by relevant published data from in vitro experiments on cancer cell viability (Ref. 10) and provided quantitative information on (i) Ncell development during incubation (48 h) and (ii) %PApoptosis and %PMitosis following cancer cell subjection to plasma-conditioned RPMI (PC-RPMI), DOX, and DOX combined with APPJ. The in silico model revealed a significant impact on Ncell, %PApoptosis, or/and %PMitosis of each “treatment case” scenario depending on dPlasma and cDOX without using any flow cytometry experiments. Thus, our simulations can help expedite experimentation and reduce related costs for exposing cancer cell death. Although our study here is focused on melanoma cell apoptosis and mitosis, other cell phenotypes (e.g., glioblastoma)23 and mechanisms (e.g., growth, migration, differentiation, etc.)25 can be specified in our model depending on the biology involved, treatment type, and application. Thus, the ABM can be expanded by adding more complexity such as RONS interaction pathways with cancer cells in different liquid media. The in silico is complementary to in vitro experiments, which is particularly useful in the plasma (cancer) medicine field for understanding in depth the pathways of cancer cell death induced by the plasma. The more evidence is provided to refine the model, and the broader the applicability of the model will be to explore more elaborate cancer treatment strategies.
See the supplementary material for a dynamic visualization of in silico Ncell development during incubation under control conditions (Animation) and graphical comparison of the in silico (our work) and in vitro results regarding Ncell, %PApoptosis, and %PMitosis values in response to PC-RPMI (Fig. S1), DOX (Fig. S2), and DOX combined with APPJs (Fig. S3). Some details on the in silico model and the APPJ electrical and optical features are also given.
This work was co-funded by the European Regional Development Fund and the Republic of Cyprus through the Research and Innovation Foundation in the framework of the project “PLASMA-TReatS-TUMORS” with Protocol No. OPPORTUNITY/0916/MSCA/0023. This project has also received funding from the European Union's Horizon 2020 research and innovation program under Grant Agreement No. 810686. We also acknowledge the financial support from the Cyprus Cancer Research Institute, as part of the project “PROTEAS” (CCRI Bridges in Research Excellence; Funding Agreement No. CCRI_2021_FA_LE_105). K. Gazeli would like thank Dr. D. Athanasopoulos for valuable discussions on the in vitro experiments.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts of interest to disclose.
Author Contributions
V.V. and G.E.G. are joint senior authors.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding authors upon reasonable request.