Anti-cancer capacity of plasma-treated PBS: effect of chemical composition on cancer cell cytotoxicity

We evaluate the anti-cancer capacity of plasma-treated PBS (pPBS), by measuring the concentrations of NO2 − and H2O2 in pPBS, treated with a plasma jet, for different values of gas flow rate, gap and plasma treatment time, as well as the effect of pPBS on cancer cell cytotoxicity, for three different glioblastoma cancer cell lines, at exactly the same plasma treatment conditions. Our experiments reveal that pPBS is cytotoxic for all conditions investigated. A small variation in gap between plasma jet and liquid surface (10 mm vs 15 mm) significantly affects the chemical composition of pPBS and its anti-cancer capacity, attributed to the occurrence of discharges onto the liquid. By correlating the effect of gap, gas flow rate and plasma treatment time on the chemical composition and anti-cancer capacity of pPBS, we may conclude that H2O2 is a more important species for the anti-cancer capacity of pPBS than NO2 −. We also used a 0D model, developed for plasma-liquid interactions, to elucidate the most important mechanisms for the generation of H2O2 and NO2 −. Finally, we found that pPBS might be more suitable for practical applications in a clinical setting than (commonly used) plasma-activated media (PAM), because of its higher stability.

(1) in which is the density of species (m -3 ) , the total number of reactions, , and , the stoichiometric coefficients at the left hand side and right hand side of the reaction and the rate of reaction (in m -3 s -1 ), given by: (2) in which is the rate coefficient (m 3 s -1 or m 6 s -1 for two-body or three-body reactions, respectively). The rate coefficients of the heavy particle reactions are either constant or dependent on the gas temperature, whereas the rate coefficients of the electron impact reactions depend on the electron temperature or the reduced electric field ⁄ (i.e., the electric field divided by the number density of all neutral species , usually expressed in Td = 10 -21 V m 2 ). The rate coefficients of the electron impact reactions are generally calculated according to the following equation: with the electron energy (usually in eV), ℎ the minimum threshold energy needed to induce the reaction, ( ) the velocity of the electrons (in m s -1 ) , ( ) the cross section of collision (in m 2 ), and ( ) the (normalized) electron energy distribution function (EEDF; in eV -1 ) calculated using a Boltzmann solver. In this work we solve the balance equations (1) of all species by means of the ZDPlaskin code, which also has a built-in Boltzmann solver, called BOLSIG+ S.1 , to calculate the EEDF and the rate coefficients of the electron impact reactions S.2 based on a set of cross sections, the plasma composition, the gas temperature and the reduced electric field (E/N). The electric field (E; in V m -1 ) is calculated from a given power density, using the so-called local field approximation S.3 : with the input power density (in W m -3 ) and the plasma conductivity (A V -1 m -1 ). The plasma conductivity is estimated at the beginning of the simulations as S.3 : with the elementary charge (1.6022x10 -19 C), , the initial electron density (in m -3 ), the electron mass (9.1094x10 -31 kg) and the collision frequency (in s -1 ) calculated using BOLSIG+. During the simulation the plasma conductivity is calculated as S.3 : = ( ) 0 (6) with the electron drift velocity (in m s -1 ), which is calculated using BOLSIG+ implemented in ZDPlaskin, and ( ) the reduced electric field at the previous time step (in V m 2 ).

Description of the plasma jet in the chemical kinetics model
In the approach of using a chemical kinetics model to simulate the kINPen plasma jet studied in this work, a cylindrical volume element is followed along the jet stream. By doing this, we assume a homogenous plasma along the radial axis (cfr. plug flow reactor). Moreover, we assume that the axial transport of mass and energy due to drift and diffusion is negligible compared to convection. Due to the very high axial flow speed (order of 10 3 cm s -1 ) compared to the radial flow speed this seems acceptable. Upon reaching the liquid substrate, the calculated gas phase densities of all plasma species are used as input for the liquid phase module. In this module, the accumulation of species in the liquid is determined by the diffusion from gas phase species into the liquid, which is based on Henry's law, as well as by the liquid-phase chemistry. This approach, which allows us to investigate the liquid chemistry using a chemical kinetics model, was introduced by Lietz et al. S.4 The general plasma jet setup, assumed in the model, is shown in Supplementary Figure S3.
Supplementary Figure S3. Plasma jet set-up used in the chemical kinetics model. The start of the simulation is 3 mm before the nozzle, which is at the tip of the inner electrode (thick gray line). The length of the visible plasma plume (indicated in purple) and the total distance between nozzle and liquid sample (both denoted as x cm) depend on the specific treatment conditions (see Table 1 in the main paper).

Gas phase module
Conceptually, a chemical kinetics model calculates the density of all species as a function of time (see equation 1). However, by assuming a certain velocity profile of the feed gas, this time can be coupled to the position of the volume element along the axis, which allows us to obtain information on the species densities as a function of distance, and thus to investigate different treatment distances, as used in the experiments. An example of the gas flow velocity profile, which decreases along the axis due to gas expansion and obstruction by the relatively stationary surrounding atmosphere, is shown in Supplementary Figure S4, for a gas flow rate of 1 slm. The initial gas flow velocity, at the nozzle, is calculated based on the flow rate of the feed gas and the dimensions of the plasma jet. Furthermore, as mentioned above, many of the gas phase reaction rate coefficients depend on the gas temperature. This means that a gas temperature profile along the plasma axis is required to calculate the exact rate of all reactions (see Supplementary Figure S4). This temperature profile is based on our experimental measurements.
Supplementary Figure S4. Plasma  Moreover, as the electron impact reactions depend on the EEDF, the reduced electric field is also required. As mentioned above, this reduced field is calculated based on the deposited power density, of which an example profile is also shown in Supplementary Figure S4. The maximum value of the power density is achieved at the tip of the powered electrode. Subsequently, the power density decreases linearly along the plasma axis, reaching zero at the end of the visible plasma plume, which is observed experimentally. This is chosen as the simulation results indicate that the densities of the excited species quickly drop to zero when the power density drops to zero, due to which the visible plasma plume would also be lost. The length of the plasma plume depends on the gas flow rate, based on our experimental observations, i.e., at 1 slm, the plasma plume propagates in general 9 mm into the surrounding atmosphere, whereas at 3 slm, the plasma plume has a length of 12 mm. In the case of 1 slm and a treatment distance of 10 mm, plasma discharges onto the liquid substrate were observed, as mentioned in the main paper. This means that under these conditions, a discharge between two electrodes occurs, (i) the electrode tip from the plasma jet and (ii) the liquid surface. Therefore, we assume the power density profile to rise again slightly upon reaching the liquid surface (i.e. at the end of the gas phase simulation). In all cases, the total deposited power equals 3.5 W, as is the case in the experimental treatments. Finally, to mimic the mixing of humid air species into the effluent of the plasma jet, these species (O2, N2 and H2O) are added into the effluent, assuming a certain air mixing rate (based on experimental data). The profiles of the ambient air species along the axis are also shown in Supplementary Figure  S4. Note that the diffusion of ambient air species only starts after 0.12 cm in the effluent. This is because it will take some time before the ambient air species are able to diffuse up to the plasma axis. The initial densities of O2, N2 and H2O inside the device (grey area in Supplementary Figure S4) originate from the impurities of the feed gas (1, 4 and 3 ppm for O2, N2 and H2O, respectively), which are taken the same as the impurities present in the feed gas used in the experimental work. The chemistry set of the gas phase reactions used in this study is largely taken from Murakami et al. S.6 However, to include additional relevant biomedically active species (e.g. H2O2, HO2, HNO3 or HNO2), we extended this chemistry set with the reactions describing the behavior of these species, adopted from the chemistry set of Van Gaens and Bogaerts S.7 , yielding a total chemistry set of 91 different species and 1390 reactions. All species included in the gas phase are shown in Supplementary Table  S1.

Supplementary
in which the first term is similar to the calculation of the gas phase species densities (conservation of mass; see above). The second term represents the diffusion of gas phase species into the liquid. In this term, Ds is the diffusion coefficient of gas phase species s, ns,g is the final gas phase density of species s and ʌ is the diffusion length of the plasma. Furthermore, fl is the fraction of the area of the plasma in contact with the liquid and Ss,l is the sticking coefficient of species s on the liquid, given by: in which hs represents the Henry constant of species s. This sticking coefficient is only used if ns,l/ns,g < hi and accounts for a diminishing rate of loss of the gas phase species into the liquid as the liquid density approaches its Henry's law equilibrium values. Finally, Vp and Vl represent the volume of the plasma and the liquid, respectively. The third term of equation 7 is only non-zero if the liquid is oversaturated (i.e. if ns,l/ns,g > hi) and represents the flux from the liquid phase into the gas phase. The Henry constants were adopted from Lietz et al. S.4 , whereas the diffusion coefficients were taken from Verlackt et al. S.5 . As mentioned before, the reaction chemistry of the liquid phase is taken from Lietz et al. S.4 and includes in total 35 species and 89 reactions. It is impossible to take into account the transportation of plasma species from the gas-liquid interface into the bulk of the liquid by means of this 0D chemical kinetics model, but in reality, the density of the short-lived reactive species, such as OH radicals, will drop quickly from the interface towards the bulk. We mimic this drop in species densities in our model by decreasing the reaction rate coefficients of the reactions involving these short-lived species with longer-lived species. This approach is based on observations from more detailed 2D fluid simulations carried out in our group, in which short-lived species react mostly at the gas-liquid interface, generating more stable species, which are then rapidly transported towards the bulk of the liquid due to convection. S.5 Finally, it is important to mention that the liquid in our model is pure water, with 4.8 ppm O2 and 8.9 ppm N2 initially dissolved into it (equilibrium values with air). The experiments were performed in a buffered solution at pH 7.3, so the concentrations of H3O + and OHin the liquid were fixed throughout the entire simulation at values which correspond to this pH.

Catalase experiments
Supplementary Figure S5