Anaïs Chanclu
email: anais.chanclu@univ-avignon.fr

Imen Ben Amor

Cédric Gendrot

Emmanuel Ferragne

Jean-François Bonastre

Automatic classification of phonation types in spontaneous speech: towards a new workflow for the characterization of speakers' voice quality

Keywords: voice quality, speaker characterization, phonation type classification, neural network, explainability

Voice quality is known to be an important factor for the characterization of a speaker's voice, both in terms of physiological features (mainly laryngeal and supralaryngeal) and of the speaker's habits (sociolinguistic factors). This paper is devoted to one of the main components of voice quality: phonation type. It proposes neural representations of speech followed by a cascade of two binary neural network-based classifiers, one dedicated to the detection of modal and nonmodal vowels, and one for the classification of nonmodal vowels into creaky and breathy types. This approach is evaluated on the spontaneous part of the PTSVOX database, following an expert manual labelling of the data by phonation type. The results of the proposed classifiers reaches on average 85 % accuracy at the framelevel and up to 95 % accuracy at the segment-level. Further research is planned to generalize the classifiers on more contexts and speakers, and thus pave the way for a new workflow aimed at characterizing phonation types.

Introduction

Voice quality is considered in the literature to have great implications for the characterization of speakers [START_REF] Gold | International practices in forensic speaker comparison[END_REF]. It can be a permanent component of a speaker's voice due to physiological particularities (mostly laryngeal and supralaryngeal) or speaker habits (sociolinguistic factors), but it may also be subject to intra-speaker variability, notably speech style or emotion [START_REF] Nolan | Forensic speaker identification and the phonetic[END_REF]. Some known examples of voice quality are nasality, tenseness/laxness, dark/clear voice, phonation type, etc.

In this study, we focus on phonation type, which is the mode of vibration of the vocal folds during voiced phonation. It is typically regarded as an articulatory continuum between closed and open glottal state, respectively for speech from creaky to modal and finally breathy voice [START_REF] Gordon | Phonation types: a cross-linguistic overview[END_REF][START_REF] Wright | Voice quality types and uses in north american english[END_REF]. Phonation type can be a phonological feature in quite a few languages. For example, Jalapa Mazatec has creaky/modal/breathy /i/, /a/, /ae/, /o/, /u/, but other languages can be mentioned such as Hindi or Tsonga. Phonation type has also been demonstrated to be used for phonetic purposes: for example, creak can be an allophone of intervocalic unvoiced stops in British English [START_REF] Garellek | Acoustic differences between english /t/ glottalization and phrasal creak[END_REF], and breathiness can be used to signal prosodic finality in French [START_REF] Smith | Marking the boundary: utterance-final prosody in french questions and statements[END_REF], hence a potential role of phonation type as speaker stylistic marker. Obviously, phonation type can be used in the phonetic characterization of a speaker or a group of speakers: women are said to exhibit a more breathy voice than men in British English [START_REF] Henton | Language, speech, and mind: Studies in honour of Victoria A[END_REF], high prestige speakers in Edinburgh have been reported to use more creaky voice [START_REF] Esling | The identification of features of voice quality in social groups[END_REF], see also [START_REF] San Segundo | The use of the vocal profile analysis for speaker characterization: Methodological proposals[END_REF] for a review.

Voice quality, including phonation type, is often studied perceptually due to limitations in phoneticians' understanding of which acoustic parameters are most relevant, and how exactly they should be measured [START_REF] Keating | Acoustic analysis of creaky voice[END_REF][START_REF] Katharina Klug | Analysing breathy voice in forensic speaker comparison. using acoustics to confirm perception[END_REF]. We aim in this article at developing a neural-network system allowing for the automatic detection of these three phonation types (modal, creaky, breathy) as they are known to be useful for the characterization of speakers' voices. The system is designed to help phoneticians learn more about phonation types. This knowledge can be easily fed back into the system, to build a virtuous circle.

Section 2 describes the system. The corpus, including the manual annotations dedicated to the targeted task, is described in Section 3. Then, Section 4 presents the experimental protocol and Section 5 the corresponding experimental results. Finally, Section 6 proposes some takeaways and future work.

PASE-MLP phonation type detection system

Our goal is to provide an automatic system able to predict the type of phonation for a given speech frame or segment. The system must be bootstrapped on a small set of monitored data because manual annotation of this type of phenomenon takes time. Its ability to use large amounts of unlabelled data later is also of great importance to improve generalization.

At the feature extraction level, the system should allow great flexibility since discovering the relevant features is a part of the task. In addition, the knowledge gathered by the expert thanks to the system should be easily reinjected into the system. The problem-agnostic speech encoder (PASE) [START_REF] Pascual | Learning Problem-Agnostic Speech Representations from Multiple Self-Supervised Tasks[END_REF] looks promising for this purpose. PASE is a variant of autoencoder [START_REF] Ng | Sparse autoencoder[END_REF][START_REF] Lu | Speech enhancement based on deep denoising autoencoder[END_REF], a kind of neural network trained without supervision to reproduce input to output. During the operating phase, a vector, the "embedding", is extracted from the weights of the central layer for each audio signal frame. PASE is trained for multitasking: the embedding extractor is optimized to achieve different objectives, materialized by a set of "workers", such as retrieving the waveform, cepstral coefficients, the phoneme or certain prosodic parameters. So, the embedding is a compressed information capable of satisfying each worker. In this work we use the PASE+ [START_REF] Ravanelli | Multi-task self-supervised learning for Robust Speech Recognition[END_REF] encoder with a set of workers we consider most useful for the task (MFCC, waveform, prosody, LIM and GIM). By selecting other workers or by defining new workers based on new knowledge, it will be possible to modify the extraction characteristics to take new knowledge into account. PASE is trained with 6609 audio files from PTSVOX [START_REF] Chanclu | Ptsvox: une base de données pour la comparaison de voix dans le cadre judiciaire," in 6e conférence conjointe Journées d'Études sur la Parole[END_REF]. The outputted frames (the embeddings) have a size of 256 coefficients, with a frame rate of 1/100s.

For the classification part of the system, we built a cascade of two binary neural network-based classifiers. The first one aims to detect modal and nonmodal frames (nonmodal frames are frames labelled creaky or breathy) in a vowel segment.

Corpus

This work is using audio data extracted from the PTSVOX [START_REF] Chanclu | Ptsvox: une base de données pour la comparaison de voix dans le cadre judiciaire," in 6e conférence conjointe Journées d'Études sur la Parole[END_REF] corpus. A specific expert-based annotation process was applied to enrich the meta information of the corpus in terms of phonation type, for the prepausal vowels. The remainder of this section shows the details of the data and the annotation process.

PTSVOX

The PTSVOX database is composed of French audio recordings. It is specifically designed to study the factors of inter-and intra-speaker variability in forensic voice comparison. These factors include speaking style (reading or spontaneous speech), recording equipment (microphone or telephone), gender, and various information about the speaker (smoking, health issues, etc.). All recordings are re-sampled to 16000Hz.

PTSVOX contains 369 speakers split in two sets: Intra composed of 24 speakers recorded multiple times and Inter, with 345 speakers recorded only once. For this work, we selected the microphone recordings of spontaneous speech from the Intra set. This set is composed of 12 female and 12 male speakers recorded once a month, using a microphone and a telephone, over a 3-month period (approximately 8 to 10 minutes per speaker whose excerpts were transcribed and segmented manually).

Voice quality annotation on prepausal vowels

More specifically, we focused on prepausal vowels because they contain linguistic key points since they occur at prosodic prominences, and they are perceptually important in communication. Prepausal vowels were perceptually annotated by one of the authors in three categories: creaky, modal, and breathy. For this annotation process, we automatically isolated vowels immediately followed by a pause for each speaker with a threshold of 50 ms for minimum vowel and pause duration. Then, by visualizing spectrograms and listening to the vowels, labels were assigned to them (see Figure 1). A single vowel could have two labels if it contained a change in voice quality. The 'other' label was assigned to ambiguous or noisy cases, which were then discarded. As there can be more than one label within a vowel, we use the term segment to refer to the labels within the vowel. The resulting corpus contains 10,361 segments for 8,889 vowels.

Experimental protocol

We decided to evaluate the two binary classifiers separately: modality detection between modal and nonmodal labels and voice quality detection between creaky and breathy labels. In the rest of this section, we present the reference system we defined, the split of the data into training and testing sets and the decision rules.

MFCC-SVM reference system

In order to serve as a baseline, we built a classical MFCC+SVM recognizer. Acoustic features are composed of 30 MFCC parameters with a 0-8kHz bandwidth, extracted using torchaudio. The window length is 20 ms with an overlap of 10 ms between 2 consecutive windows, giving a frame rate of 100 frames per second. We used a support vector machine (SVM) with a radial basis function. Similarly to our MLP system, the models were trained on a balanced and randomized amount of data between classes.

Training and test subsets

For each gender, we randomly selected two speakers and used the corresponding audio data to form a fully separated test set. The selected speakers for the test set are LG001 and LG009 for female speakers and LG012 and LG024 for male speakers. 20 % of the remaining data are then randomly selected to form a validation set, and the rest is dedicated to the training set. The speakers of the test set are never seen in the training or in the validation sets. The train and validation sets contain recordings spoken by the same speakers. 1 presents the number of annotated prepausal vowel segments and the corresponding duration (in number of frames) for both sets and vowel label.

Scoring rules

The predictions are done frame by frame for both tasks (modal/nonmodal and creaky/breathy) and both recognizers, with a rate of 100 frames per second (recall that the parameterization steps of the two recognizers used a signal window longer than 1/100s to compute a given frame). Decisions were made using an simple rule: if the score was higher than 0.5, the rec-ognizer answered yes, and no otherwise (there was no threshold setting at all). Performance at the segment-level was also computed. To predict voice quality classification at the segmentlevel, we first performed a prediction at the frame-level. Then we computed the mean prediction over all frames of the segment before applying the same decision rule as previously.

Results

We first present the results for the modal/nonmodal task before moving on to those of the creaky/breathy classification task. Performance is presented using confusion matrices for a clearer visualization of class imbalance in the test set. Performance on the test set only (i.e. using data coming from speakers never seen during the training) is shown.

Modal/nonmodal classification

Table 2 presents the performance of the modal/nonmodal classifiers at the frame-level. 87 % of the modal frames and 76 % of the nonmodal frames are well classified with the PASE-MLP system (κ = 0.62). The reference system shows lower performance with, respectively, 79 % and 72 % (κ = 0.49). Table 3 shows the performance for the same task but gathered at the segment-level. Overall, the same difference between our system PASE-MLP and the baseline appears. As expected, the general level of performance is significantly higher than for the framelevel, with up to 93 % of good classification for modal frames (κ = 0.76). We can notice the performance of the baseline system at the segment-level is close to the performance of our PASE-MLP system at the frame-level (κ = 0.65).

Creaky/breathy classification

Table 4 presents the performance of the creaky/breathy classifiers at the frame-level. The baseline system has a good recognition rate for breathy frames (88 %) but shows a clear weakness for creaky frames with 67 % of correct recognition (κ = 0.51).

Our PASE-MLP approach outperforms the baseline with 85 % for creaky frames and 88 % for breathy frames (κ = 0.73). This favours the cascade architecture: it seems easier to separate breathy and creaky frames than modal and nonmodal frames, although additional experiments are required before a definitive conclusion can be reached. Table 5 shows the results of similar experiments to Table 4 but at the segment-level. The performance confirms our previous findings. The performance at the segment-level for this task with our PASE-MLP system appears very promising, with about 95 % of correct classification (κ = 0.88, and κ = 0.67 for MFCC-SVM).

Discussion and conclusion

This article aims to lay the foundations of a new workflow dedicated to voice quality, and more precisely to phonation type. The final objective is to propose an automatic tool to help phoneticians to improve their understanding of phonation phenomena.

We started here at a low level by defining a new neural network solution to automatically detect three phonation types (modal, creaky and breathy) on prepausal vowels. We implemented a cascade of two binary classifiers, based on a representation learning approach, PASE, for feature extraction and a classical MLP for classification. These choices were driven by two wishes: to be able to generalize the knowledge gathered from the -limited in size-manual annotation provided by phoneticians, and to embed easily the phoneticians understanding in the system, following a virtuous circle model. The proposed solution yielded very satisfactory performance, with on average ≈84 % of correctly labelled frames for both classifiers. The performance reached ≈91 % when the decision was taken at the segment level. Our neural network-based solution, PASE-MLP, clearly outperformed a classical MFCC-SVM approach (≈75 % on average at the frame level for the latter, to be compared to ≈85 % for the PASE-MLP system). More importantly, the performance of PASE-MLP system appeared more stable than the MFCC-SVM solution.

The results open the door for further investigation. First, the PASE-MLP system should be generalized to more phonetic contexts and speakers. Our next step is to use transfer learning [START_REF] Pan | A survey on transfer learning[END_REF] to pretrain both PASE and MLP modules (it will use the current system to propose initial labels as suggested in [START_REF] Gresse | Learning Voice Representation Using Knowledge Distillation For Automatic Voice Casting[END_REF]). It will allow us to apply the system on all the PTSVOX data as well as on other, larger, databases. Second, the system will provide expert phoneticians with the best subset of examples for perceptual study and the possibility to incorporate their feedback. Finally, we wish to help expert phoneticians to study speaker phonotypes on a large scale, including intra-speaker variability, thanks to the automatic tools we will provide.

Figure 1 :

 1 Figure 1: Labelling example of a creaky vowel

Table 1 :

 1 Number of segments and duration (in 1/100s frames) for the train and test subsets

		Train & Val		Test
		Seg.	Frames Seg. Frames
	modal	4,090 89,478	687	15,773
	nonmodal 2,220 33,554	464	7,410
	creaky	1,088 17,441	251	3,569
	breathy	1,132 16,113	213	3,841
	Table

Table 2 :

 2 Confusion matrices between modal and nonmodal frames in the test set for both systems

			PASE-MLP
		modal	nonmodal	Total
	modal	13,640 (87 %) 2,133 (13 %) 15,773
	nonmodal	1,768 (24 %)	5,642 (76 %)	7,410
			MFCC-SVM
	modal	12,430 (79 %) 3,343 (21 %) 15,773
	nonmodal	2,070 (28 %)	5,334 (72 %)	7,410

Table 3 :

 3 Confusion matrices between modal and nonmodal segments in the test set for both systems

		PASE-MLP
		modal	nonmodal	Total
	modal	638 (93 %)	49 (7 %)	687
	nonmodal	80 (17 %)	384 (83 %)	464
		MFCC-SVM
	modal	595 (87 %)	92 (13 %)	687
	nonmodal 102 (22 %) 362 (78 %)	464

Table 4 :

 4 Confusion matrices between creaky and breathy frames on the test set for both systems

			PASE-MLP
		creaky	breathy	Total
	creaky	3,048 (85 %)	521 (15 %)	3,569
	breathy	465 (12 %)	3,376 (88 %) 3,841
		MFCC-SVM
	creaky	2,376 (67 %) 1,193 (33 %) 3,569
	breathy	586 (15 %)	3,255 (85 %) 3,841

Table 5 :

 5 Confusion matrices between creaky and breathy segments on the test set for both systems

		PASE-MLP
		creaky	breathy	Total
	creaky	235 (94 %)	16 (6 %)	251
	breathy	11 (5 %)	202 (95 %)	213
		MFCC-SVM
	creaky	193 (77 %)	58 (23 %)	251
	breathy	19 (9 %)	194 (91 %)	213

Acknowledgements

The research reported here was supported by the ANR-17-CE39-0016 VoxCrim project.


J Ramos-Méndez

J Laverne

N Domínguez-Kondo

J Milligan

V Štěpán

Laverne Ja

Domínguez-Kondo N

W-G Shin

TOPAS-nBio validation for simulating water radiolysis and DNA damage under low-LET irradiation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

simulations but disagreed by 12% with published direct measurements. Improvement of the spatial resolution of the DNA damage model might mitigate such disagreement.

In conclusion, with these improvements, Geant4-DNA/TOPAS-nBio provides a fast, accurate, and user-friendly tool for simulating DNA damage under low LET irradiation.

Introduction.

A reliable way to study the underlying channels of radiobiological damage is mechanistic modeling. The modeling task is accomplished with Monte Carlo track-structure (MCTS) codes that have the capability of combining complex DNA geometry models with the stochastic processes of the interaction of ionizing radiation with matter and the subsequent non-homogeneous chemistry processes to produce initial DNA damage. MCTS codes specific for radiobiological applications include KURBUC [START_REF] Nikjoo | Radiation track, DNA damage and response-a review Reports[END_REF], PARTRAC [START_REF] Dingfelder | Electron inelastic-scattering cross sections in liquid water[END_REF], Geant4-DNA (Incerti et al 2010a, 2010b[START_REF] Bernal | Track structure modeling in liquid water: A review of the Geant4-DNA very low energy extension of the Geant4 Monte Carlo simulation toolkit[END_REF], RITRACKS [START_REF] Plante | Considerations for the independent reaction times and step-bystep methods for radiation chemistry simulations[END_REF], RADAMOL [START_REF] Štěpán | RADAMOL tool: Role of radiation quality and charge transfer in damage distribution along DNA oligomer Guest editors[END_REF], and TOPAS-nBio [START_REF] Schuemann | TOPAS-nBio: An Extension to the TOPAS Simulation Toolkit for Cellular and Sub-cellular[END_REF]. The reported results obtained with such codes demonstrate the usability of the Monte Carlo method to quantify DNA damage-induced at early stages following irradiation from first principles.

MCTS codes, being flexible mechanistic tools initially developed for liquid water calculations, have demonstrated acceptable accuracy in estimating DNA damage induced by ionizing radiation from first principles. However, many radiobiology-oriented MCTS codes overestimate the yields of • OH radical produced by fast electrons in liquid water within the nanosecond time scale [START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF][START_REF] Uehara | Monte Carlo simulation of water radiolysis for low-energy charged particles[END_REF], Ramos-Méndez et al 2018[START_REF] Boscolo | TRAX-CHEM: A pre-chemical and chemical stage extension of the particle track structure code TRAX in water targets[END_REF]. This inaccuracy might potentially mislead the interpretation of calculated DNA damage yields and hence requires benchmarking. In this work, we aim to reconcile MCTS by comparing simulated with experimentally measured yields from the literature for radiolytic species produced by fast electrons at the earliest times.

Experimental validation is paramount to determine the accuracy of the Monte Carlo method. A common approach among MCTS codes assumes that, from the physics perspective, liquid water is sufficient to represent biological tissue (Friedland et al 2017[START_REF] Nikjoo | Radiation track, DNA damage and response-a review Reports[END_REF], Sakata et al 2019, Zhu et al 2020b). For this reason, MCTS validation relies upon the comparison between calculated quantities with measurements performed in gas or liquid water (see e.g., [START_REF] Burigo | Distributions of deposited energy and ionization clusters around ion tracks studied with Geant4 toolkit[END_REF][START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF][START_REF] Pimblott | Stochastic Simulation of the Electron Radiolysis of Water and Aqueous Solutions[END_REF]). Under low linear energy transfer (LET) irradiation, e.g., 0.4 keV/µm track-averaged LET for 137 Cs [START_REF] Bruce W R, Pearson | The Linear Energy Transfer Distributions Resulting from Primary and Scattered X-Rays and Gamma Rays with Primary HVL's from 1.25 mm Cu to 11 mm[END_REF], the highest contribution to DNA lethal damage is caused by reactions with hydroxyl radicals produced in the radiolysis of water (e.g. > 70% for 60 Co, estimated with DNA plasmids pBR322 at scavenging capacities below ~5 x 10 8 s -1 ) [START_REF] Klimczak | Irradiation of plasmid and phage DNA in water-alcohol mixtures: Strand breaks and lethal damage as a function of scavenger concentration Int[END_REF]. Therefore, thorough validation of the models used by MCTS to simulate the chemical stage in the radiolysis in water is crucial. This task has been accomplished typically by comparing calculated time-dependent yields in pure liquid water with measurements performed in solutes at different scavenging capacities. Subsequently, MCTS codes have been tuned to match the experimental yields at the earliest times available to date.

The accuracy of experimental measurements of radiolytic yields continuously improves. Nowadays, the most accurate picosecond-level measurements of radiolytic yields produced in liquid water by fast electrons have been performed with electron pulse radiolysis. The most recent measurements have determined reference G values (number of chemical species created or lost per 100 eV of energy deposit) for hydroxyl radicals ( • OH) of 4.8 ± 0.2 molec./100 eV and solvated electron (e - aq ) yields of 4.2±0.2 molec./100 eV after just 7 ps (El Omar et al 2011) [START_REF] Wang | Timedependent yield of the hydrated electron and the hydroxyl radical in D2O: A picosecond pulse radiolysis study[END_REF]. These values are significantly lower than previous values used to calibrate Monte Carlo simulations (see e.g., [START_REF] Tomita | Monte Carlo simulation of physicochemical processes of liquid water radiolysis[END_REF][START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF][START_REF] Pimblott | Reconciliation of Transient Absorption and Chemically Scavenged Yields of the Hydrated Electron in Radiolysis[END_REF][START_REF] Uehara | Monte Carlo simulation of water radiolysis for low-energy charged particles[END_REF]) and call for revisiting the input parameters of MCTS codes. A re-evaluation of such parameters has been mentioned for the radiation chemistry code IONLYS-IRT [START_REF] Sanguanmith | Time-dependent yield of OH radicals in the low linear energy transfer radiolysis of water between 25 and 350[END_REF][START_REF] Sultana | Yields of primary species in the low-linear energy transfer radiolysis of water in the temperature range of 25-700[END_REF]; however, specific details of parameter adjustment were not reported.

On the other hand, to verify the accuracy of MCTS for radiobiology applications, sophisticated DNA models have been developed based on the whole cellular nucleus to identify and quantify clustered DNA strand breaks sites [START_REF] Nikjoo | Radiation track, DNA damage and response-a review Reports[END_REF], Friedland et al 2017[START_REF] Štěpán | RADAMOL tool: Role of radiation quality and charge transfer in damage distribution along DNA oligomer Guest editors[END_REF][START_REF] Meylan | Simulation of early DNA damage after the irradiation of a fibroblast cell nucleus using Geant4[END_REF], Lampe et al 2018b, Sakata et al 2019, Zhu et al 2020a). The spatial resolution of experimental measurements at the cellular level is insufficient to validate Monte Carlo simulation outputs directly; for example, experimental techniques can still not resolve spatial distributions of individual DNA single-strand breaks (SSB) and double-strand breaks (DSB), and there is insufficient data about the structure and the radiation chemistry of chromatin. Thus, different assumptions have been made by different research groups to reconcile their MCTS results with experimental data, which impacts the accuracy of the codes. In particular, the scavenging capacity of the biological environment has been mimicked by limiting the time domain to a few nanoseconds. [START_REF] Moiseenko | The cellular environment in computer simulations of radiation-induced damage to dna[END_REF] showed that this approach introduces 20%-30% differences in strand break yields compared to the explicit simulation of scavenger molecules uniformly distributed around the radiation spurs. Besides, the simulation of the chemical and biological mechanisms occurring within a cell is not a trivial task [START_REF] Wardman | Radiotherapy Using High-Intensity Pulsed Radiation Beams (FLASH): A Radiation-Chemical Perspective[END_REF], and MCTS simulations have not yet utilized the corresponding models to simulate more reliable cellular conditions.

The selection of a less complex biological system than a cell would assist in providing relevant data to evaluate the accuracy of MCTS. Plasmids of DNA are a convenient example. Experiments consisting of the irradiation of plasmids can be performed in a well-controlled environment with solutes of various scavenging capacities [START_REF] Milligan | Characterization of the reaction rate coefficient of DNA with the hydroxyl radical[END_REF]. Thus, under low LET irradiation, direct comparison with plasmid experiments provides a suitable way to validate directly MCTS radiochemistry capabilities. Comparison of MCTS with plasmid experiments considering direct and indirect damages has been performed in the past [START_REF] Tomita | Monte Carlo simulation of DNA strand-break induction in supercoiled plasmid pBR322 DNA from indirect effects[END_REF] (Fulford et al 2001) [START_REF] Edel | Computer simulation of strand break yields in plasmid pBR322: DNA damage following 125I decay[END_REF].

In this work, we focus on the comparison of calculated versus published measured data to evaluate the accuracy of TOPAS-nBio. Parameters for modeling the radiation chemistry process of TOPAS-nBio were revisited, driven by the latest published experimental measurements of G values of radiation yield for fast electrons at the picosecond level. In that way, a reconciliation between measured and Monte Carlo modeled G values in liquid water was achieved. Then, we validated TOPAS-nBio for simulation of indirect damage of DNA. To this end, experimental setups of plasmids under low-LET irradiation reported in the literature were replicated with our validated Monte Carlo track-structure tool. Our benchmarked code provides a tool over which extended features regarding effects of temperature, compaction of DNA, oxygen concentration and other radiation chemistry processes may be implemented and evaluated.

2 Materials and methods.

The physics list comprises models describing the discrete transport of electrons in liquid water with no production cuts for secondary electrons, and all the interaction events are explicitly simulated. The selected models have been described extensively elsewhere (Incerti et al 2018, Shin et al 2019, Ramos-Méndez et al 2020). In brief, the physics list, encapsulated in the "G4EmDNAPhysics_option2" constructor, includes an elastic scattering model based on the partial wave theory and an inelastic scattering model based on the formalism of the complex dielectric response function of liquid water. For the latter model, four ionization shells and five discrete electronic excitation states are considered. In Geant4-DNA, after an ionization has taken place, the energy of the emitted secondary electron is calculated from the differential (in energy transfer) ionization cross section, which is described by the so-called Born ionization model. This model applies to weakly bound electrons in the liquid water molecule and adopts the dielectric response function formalism. In the case of K-shell ionization (of the oxygen atom), the atomic model Binary-encounter-Approximation-with-Exchange (BEAX) is used. Both models are described in more detail in [START_REF] Incerti | Geant4-DNA example applications for track structure simulations in liquid water: A report from the Geant4-DNA Project[END_REF] and references therein. Models to simulate vibrational excitation and electron attachment processes are also included in the constructor.

The reaction kinetics is calculated with the independent reaction times method, IRT [START_REF] Tachiya | Theory of diffusion-controlled reactions: Formulation of the bulk reaction rate in terms of the pair probability[END_REF][START_REF] Clifford | Stochastic Models of Multi-species Kinetics in Radiation-induced Spurs[END_REF], Green et al 1990[START_REF] Pimblott | Stochastic models of spur kinetics in water[END_REF]. The implementation in TOPAS-nBio is described elsewhere [START_REF] Schuemann | TOPAS-nBio: An Extension to the TOPAS Simulation Toolkit for Cellular and Sub-cellular[END_REF], Ramos-Méndez et al 2020). The simulation of scavengers is performed with the continuum approximation [START_REF] Pimblott | Stochastic models of spur kinetics in water[END_REF] where it is assumed that the scavenging molecules are uniformly distributed in the background. Then, the probability of chemical species potentially reacting with the background at time t is described by an exponential distribution given by 1 -exp(-k[B] t), where the product of reaction rate k and scavenger concentration [B], k[B], is the scavenging capacity of the background. Instantaneous scavenging within reaction distance R, at the time of chemical species creation, is performed with a probability given by exp(-4πR 3 [B]/2) [START_REF] Pimblott | Stochastic models of spur kinetics in water[END_REF]. The reactions and rate constants used in this work are presented in table 1, obtained from [START_REF] Pimblott | Investigation of various factors influencing the effect of scavengers on the radiation chemistry following the high-energy electron radiolysis of water[END_REF]). These values originated from the National Institute of Standards and Technology (NIST) database [START_REF] Buxton | Critical Review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radicals in Aqueous Solution[END_REF], where each rate constant was evaluated for accuracy and consistency.

Validation of TOPAS-nBio for water radiolysis simulations for fast electrons.

To validate TOPAS-nBio, we revisit Geant4-DNA (version 10.6.p03) parameters for the simulation of the pre-chemical stage of water radiolysis previously reported [START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF][START_REF] Karamitros | Modeling Radiation Chemistry in the Geant4 Toolkit[END_REF], Ramos-Méndez et al 2018, Shin et al 2019). Prior to the dissociation of the ionized water molecules (H 2 O + ), hole migration by electron correlation and electron relaxation (charge migration) results in a displacement of these molecules and their dissociation products from the place of energy transfer [START_REF] Ogura | Positive hole migration in pulse-irradiated water and heavy water[END_REF][START_REF] Despré | Attosecond hole migration in benzene molecules surviving nuclear motion[END_REF][START_REF] Kuleff | Core Ionization Initiates Subfemtosecond Charge Migration in the Valence Shell of Molecules[END_REF]. In Geant4-DNA, this process is considered by adding an isotropic displacement from where the ionization event occurred. The distance is sampled from a normal distribution with a given rootmean-square deviation (RMS). Due to the lack of measured data at the pre-chemical stage, the root-mean square (RMS) values (inherited from the PARTRAC code) are set to match measured yields of • OH radical available in 2009 [START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF]. In Table 2, the distances used in Geant4-DNA are described. In this work, the selection of empirical values and is driven 𝑅𝑀𝑆 𝐻 2 𝑂 + 𝑅𝑀𝑆 𝐻 2 𝑂 * by experimental G values for • OH radicals obtained at the picosecond time range for fast electrons reported in [START_REF] Wang | Timedependent yield of the hydrated electron and the hydroxyl radical in D2O: A picosecond pulse radiolysis study[END_REF]. We increased the displacement distance in multiples of the mean separation between water molecules, 0.31 nm [START_REF] Perkins | Protein volumes and hydration effects. The calculations of partial specific volumes, neutron scattering matchpoints and 280-nm absorption coefficients for proteins and glycoproteins from amino acid sequences[END_REF], to minimize the discrepancy with measured data. The branching ratios are kept unchanged, as reported in table 4 in [START_REF] Shin | Evaluation of the influence of physical and chemical parameters on water radiolysis simulations under MeV electron irradiation using Geant4-DNA[END_REF] and originally reported in [START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF].

The thermalization of sub-excited electrons is simulated using a "one-step model" with displacement parameters obtained from the literature. The transport of each individual subexcited electron is stopped once its kinetic energy is reduced to <10 eV. Later, in this one-step model, a solvated electron is placed at a distance (around the last interaction point) randomly sampled from a Gaussian distribution, with a standard deviation given by the data reported in [START_REF] Ritchie | Interactions of Low-Energy Electrons with Condensed Matter: Relevance for Track Structure Computational Approaches in Molecular Radiation Biology[END_REF], with the correction factor described in [START_REF] Shin | Evaluation of the influence of physical and chemical parameters on water radiolysis simulations under MeV electron irradiation using Geant4-DNA[END_REF]. The data from [START_REF] Ritchie | Interactions of Low-Energy Electrons with Condensed Matter: Relevance for Track Structure Computational Approaches in Molecular Radiation Biology[END_REF], provided the best agreement with measured data of the time-dependent G-values for solvated electrons from (El Omar et al 2011) when using the physics list selected in this work. This model and other approaches to simulate the thermalization of solvated electrons have been discussed in detail elsewhere [START_REF] Shin | Evaluation of the influence of physical and chemical parameters on water radiolysis simulations under MeV electron irradiation using Geant4-DNA[END_REF].

Table 1 List of reaction rate constants used for the simulation of G values obtained from [START_REF] Buxton | Critical Review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radicals in Aqueous Solution[END_REF] [START_REF] Milligan | Characterization of the reaction rate coefficient of DNA with the hydroxyl radical[END_REF] (Pastina and LaVerne 1999) The simulation setup consists of a homogenous liquid water box (1 g/cm 3 density) irradiated with monoenergetic electrons of 1 MeV. This setup is commonly used to simulate MCTS codes for fast electrons, as reported elsewhere [START_REF] Pimblott | Stochastic Simulation of the Electron Radiolysis of Water and Aqueous Solutions[END_REF][START_REF] Uehara | Monte Carlo simulation of water radiolysis for low-energy charged particles[END_REF], Ramos-Méndez et al 2018). Specifically, an isotropic electron source is positioned in the center of a cubic water phantom of 1 cm side. When the primary electron has lost more than 10 keV, the tracking of the electron is stopped, and it is removed from the simulation. The secondary electrons are then followed until their kinetic energy is reduced at or below 10 eV, the low energy limit of the ELSEPA model [START_REF] Shin | Development of a new Geant4-DNA electron elastic scattering model for liquid-phase water using the ELSEPA code[END_REF]. After that, thermalization is simulated with the one-step model. The total energy deposited in the aqueous medium thus corresponds to the energy lost by the primary electron, which is at least 10 keV [START_REF] Karamitros | Extension de l'outil Monte Carlo généraliste Geant4 pour la simulation de la radiolyse de l'eau dans le cadre du projet Geant4-DNA[END_REF]. Immediately, the pre-chemical stage takes place, and the initial position of radiolytic products are obtained and input to the IRT method for the simulation of the chemical stage up to 10 µs. (1) If e - aq are produced after an autoionization event, then its position is sampled using the onestep thermalization model. See the text.

(2) Either value is selected randomly with 50% probability.

(3) Positioned in the opposite direction from the first • OH A comparison with reference data for H 2 , H •, and H 2 O 2 was performed. For H 2 , data was obtained from (Pastina et al 1999). The authors stated that reported G values are suitable for direct comparison with Monte Carlo simulations without scavengers. For H 2 O 2 and H • , the scavenger system used in the experiments was replicated in our simulations. G values of H 2 O 2 measured in solutions containing NO 3 -at a concentration of 25 mM, and CH 3 OH at concentrations ranging from 10 -3 -10 M, were obtained from (Hiroki et al 2002). For H • , measured G values in solutions containing 1 mM of Br, 1 mM of NO 3 -and HCO 2 -at a concentration ranging from 10 -2 -1 M were obtained from (Huerta Parajon et al 2008). For the latter set of data, replicating the experiment, the calculated yield of H • was estimated from the total yield of H 2 calculated with HCO 2 -subtracted from the total yield of H 2 calculated without HCO 2 -(Huerta Parajón 2010). The reactions used for simulations in the presence of scavengers are shown in table 1. The irradiation setup was the same as that used for pure liquid water simulations.

2.2 Validation of TOPAS-nBio for damage in supercoiled DNA plasmid.

The validation of DNA damage under low-LET irradiation was performed for the experimental conditions from plasmid DNA in aerobic aqueous solution irradiated with 137 Cs γrays, as reported in [START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF][START_REF] Milligan | Yield of single-strand breaks due to attack on DNA by Scavengerderived radicals[END_REF]. A two-stage simulation was performed using condensed-history and track-structure Monte Carlo with TOPAS and TOPAS-nBio, respectively, as described below.

The first stage was used to determine the energy spectrum, at electronic equilibrium conditions, of secondary electrons set in motion by 137 Cs γ-rays interacting in a water phantom. For that, the TOPAS tool [START_REF] Perl | TOPAS: An innovative proton Monte Carlo platform for research and clinical applications[END_REF], Faddegon et al 2020) version 3.5 was used to simulate a concentric system of two homogenous spherical water phantoms of 5 cm and 10 cm radius, respectively. The physics list used the constructor "G4EmStandardPhysics_option4" with a production cut for secondary electrons of 0.05 mm, extensively benchmarked for applications in medical physics [START_REF] Arce | Report on G4-Med, a Geant4 benchmarking system for medical physics applications developed by the Geant4 Medical Simulation Benchmarking Group[END_REF]. An isotropic point source of monoenergetic γ-rays (662 keV) was positioned at the center, as shown in figure 1. We obtained the vertex kinetic energy spectrum (at the position of creation) of those secondary electrons that reached a phase space tallied on the surface of a 5 cm radius sphere. In total, 4 x 10 8 γ-rays were simulated.

The second stage consists of track-structure Monte Carlo simulations with TOPAS-nBio. The following assumptions are made for simulating the interactions of secondary electrons and reactions of chemical species with DNA molecules. All the medium, including the region occupied by the DNA geometry is made of liquid water. For biological targets, a canonical double helix B-DNA configuration is considered. The DNA sugar-phosphate groups and nucleobases, represented as cut spheres in figure 1, are the main target of the radiation or chemical species to cause DNA damage. The spatial coordinates of the centroids of the semi-spheres corresponding to sugar-phosphate groups are included in the IRT method for the realization of reactions with the chemical species produced in the water radiolysis process. The coordinates are determined by the supercoiled path, which defines the plasmid DNA model (see below). DNA nucleobases are not included in the reaction kinetics because observed rate constants between chemical species and DNA are used [START_REF] Tomita | Monte Carlo simulation of DNA strand-break induction in supercoiled plasmid pBR322 DNA from indirect effects[END_REF][START_REF] Perry C C, Ramos-Méndez | DNA condensation with a boron-containing cationic peptide for modeling boron neutron capture therapy[END_REF], but they are present for ionization/excitation interactions (see below). For the simulated time domain of the chemical stage (10 -12 -10 -4 s), the plasmid DNA model is assumed static. Accumulation of energy deposition of at least 17.5 eV in the sugar-phosphate volumes is considered to register an SSB from direct physical interactions (see (Lampe et al 2018a, Zhu et al 2020b) and references therein). Chemical species originated from radiolysis occurring within the regions occupied by the DNA semi-spheres are not included in the reaction kinetics and are eliminated upon creation, i.e., no DNA radiolysis is simulated.

For modeling indirect damage, a SSB is registered with specific efficiency after the reaction between • OH and DNA (as a whole) occurred [START_REF] Önal | γradiolysis of poly(a) in aqueous solution: Efficiency of Strand break formation by primary water radicals Int[END_REF][START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF][START_REF] Klimczak | Irradiation of plasmid and phage DNA in water-alcohol mixtures: Strand breaks and lethal damage as a function of scavenger concentration Int[END_REF]. The reaction rate coefficient for the •OH+DNA reaction (table 1) depends on the scavenging capacity of the irradiated environment and is obtained from measured data reported elsewhere [START_REF] Milligan | Characterization of the reaction rate coefficient of DNA with the hydroxyl radical[END_REF]. For H • and e - aq , constant reaction rates are used as obtained from [START_REF] Buxton | Critical Review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radicals in Aqueous Solution[END_REF]. We estimate the DNA strand break efficiencies of • OH and H • by minimization. We compare calculated to measured data applying a minimization algorithm using the Nelder-Mead method [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] to find the best efficiency values. For e - aq , no strand breaking in DNA is counted as it has not been observed experimentally [START_REF] Jones | Kinetics of radiation-induced Strand break formation in singlestranded pyrimidine polynucleotides in the presence and absence of oxygen; a timeresolved light-scattering study[END_REF], but reactions with DNA are included as e - aq binds efficiently with nucleobases [START_REF] Kumar | Reaction of electrons with dna: Radiation damage to radiosensitization[END_REF].

The DNA model is wrapped around supercoiled paths to reconstruct pUC18 plasmids (length of 2686 base-pairs length) utilizing DNAfabric [START_REF] Meylan | Geant4-DNA simulations using complex DNA geometries generated by the DnaFabric tool[END_REF]. In a separate simulation, the geometric model of the pUC18 plasmid (Yanisch-Perron et al 1985) is constructed by the worm-like chain method and subsequent smoothing to base-pair level. First, the supercoiling of a circular plasmid represented by 91 linear segments is modeled using the elastic worm-like chain code developed by the group of A. Vologodskii [START_REF] Vologodskii | Conformational and thermodynamic properties of supercoiled DNA[END_REF]Cozzarelli 1994, Huang et al 2001). The temperature is set to 298 K, the superhelical density to -0.06, and the step length to 500 ps. Resulting coarse configurations are converted to sets of 2686 equidistant points along a smooth path using an in-house code implementing the approach of Kümmerle and Pomplun [START_REF] Kümmerle | A computer-generated supercoiled model of the pUC19 plasmid[END_REF]. One sample supercoiled plasmid configuration is selected for the Monte Carlo simulation of DNA damage. The plasmid is positioned multiple times inside a spherical water phantom of 0.5 µm radius. The positions and orientations of the plasmids are uniformly random. The sphere is centered in a cubic water box of 2 µm side. A volumetric electron source is uniformly distributed in the cubic box, including the sphere. The initial electron spectrum is calculated in the first stage using a condensed-history Monte Carlo, where the initial directions are set to be isotropic (figure 1). In the experiment [START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF], the DNA is dissolved in an aerated solution containing DMSO. Thus, we simulated the scavenging behavior of that solution using the reaction rates shown in Table 1 for e - aq and H • with O 2 , using a concentration of 21% O 2 (0.27 x 10 -3 mol dm -3 ). We report the G value of SSBs using the same units of the experimental data (µmol J -1 ) as a function of scavenging capacity for DMSO. In this way, measured results from (Perry et al 2021) performed with different scavenger can be included. The DMSO concentrations range from 0.5 x 10 -4 to 1 mol dm -3 and the DNA concentration is 50 µg mL -1 . We estimated that nine pUC18 plasmids in the spherical phantom of 0.5 µm radius are equivalent to approximately 50 µg mL -1 of DNA concentration (50 x 10 -6 g mL -1 x 10 3 L -1 x mL /(2686 bp x 650 g mol -1 bp -1 )] x N av x 5.24 x 10 -16 L 9). Results for calculated SSB yields as a function of DNA concentration (from 10 -200 ≈ µg mL -1 ) for a concentration of 10 -3 mol dm -3 DMSO are also reported. Finally, the yield of SSB and DSB as a function of • OH scavenging capacity (from 7.1 x 10 5 -7.1 x 10 9 s -1 ) are reported. A DSB was scored when there was no more than a 10 base-pair distance between at least two SSB's occurring in opposite DNA strands. The total number of simulation jobs with independent random number seeds were 200 to 500, achieving statistical uncertainties from 0.7% to 5% (one standard deviation) for the simulation with the lowest and highest DMSO concentration, respectively. On average, 1150±80 primary histories are generated in each simulation job to achieve an absorbed dose of 30 Gy in the spherical phantom. In a systematic study (not shown) we found that 30 Gy absorbed dose provided a reasonable computing memory and speed, resulting from the number of primary yields used in the IRT and their processing time. The simulations take between 3 to 6 minutes to complete on a single core of a 2.7 GHz 12-Core Intel Xeon E4. 3 Results.

3.1 Time-dependent G values for fast electrons.

It was found that a distance of four water molecule's mean separation distances, 1.24 nm, for both and reproduced the measured time-dependent G values for • OH 𝑅𝑀𝑆 𝐻 2 𝑂 + 𝑅𝑀𝑆 𝐻 2 𝑂 * radicals. The G values using that value are shown in figure 2. In the top row of the figure, the calculated G value in pure liquid water for • OH, e - aq and H 2 are displayed in individual panels, along with experimental data. The agreement for the three sets of data along all the time domains was within one standard deviation of experimental measurements. In the bottom row of figure 2, the G value for H 2 O 2 and H • is displayed in each panel. As depicted, the time-dependent G values calculated in pure liquid water (solid line) did not reproduce the behavior of the experimental yields for H 2 O 2 (empty circles) and H • (filled triangles). However, the simulations considering the scavengers used in the experiments (squares connected by dashed lines) agreed within one standard deviation of the experimental data. (solid lines). Measured data for pUC18 is from ◼ [START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF], ◆ [START_REF] Milligan | Yield of single-strand breaks due to attack on DNA by Scavengerderived radicals[END_REF], and  [START_REF] Milligan | Characterization of the reaction rate coefficient of DNA with the hydroxyl radical[END_REF]. Measured data for pUC19 ▲ is from [START_REF] Perry C C, Ramos-Méndez | Boronated Condensed DNA as a Heterochromatic Radiation Target Model Biomacromolecules acs[END_REF]. Percentage differences (☐) between calculated to measured data are shown in the bottom figure. The dotted lines are margins from experimental uncertainty. Error bars are displayed when bigger than the symbol, represent statistical uncertainties from Monte Carlo simulations, one standard deviation.

3.2 DNA damage in plasmid geometries.

In figure 3, calculated and measured SSB yields as a function of • OH scavenging capacity are shown. The estimated efficiencies for • OH and H • are 24% and 0.5%, respectively. These values are used in all the following results. As depicted, the calculated SSB yields reproduced the measured data from pUC18 irradiations along with the considered DMSO concentrations reasonably well. The yields included the G(SSB) produced by direct effects, which resulted in 1.77±0.01 x 10 -6 µmol/J for the geometrical DNA model used in this work. Figure 4 shows G(SSB) as a function of the DNA concentration. The slopes of linear regression fitting to both calculated and measured data agreed within 1%±0.8%. Finally, calculated DBS as a function of • OH scavenging capacity reproduced the behavior of measured yields well as shown in figure 5. [START_REF] Klimczak | Irradiation of plasmid and phage DNA in water-alcohol mixtures: Strand breaks and lethal damage as a function of scavenger concentration Int[END_REF], pBR322 [START_REF] Tomita | Strand Break Formation in Plasmid DNA Irradiated in Aqueous Solution: Effect of Medium Temperature and Hydroxyl Radical Scavenger Concentration[END_REF]. ∇ 4 Discussion.

In this work, TOPAS-nBio was applied to simulate water radiolysis and DNA damage under low-LET irradiation. The simulation of the reaction kinetics was performed with the IRT method, providing a fast and reliable tool to assist in investigating the biological effect of the interaction of ionizing radiation at the early stages.

Calculated G values at the picosecond stage agreed within experimental errors of published direct measurements of • OH and e - aq performed in pure liquid water. The temporal evolution of the G values simulated up to the microsecond time stage was also well reproduced by TOPAS-nBio for • OH, e - aq, and H 2 . To obtain such an agreement we adjusted the RMS of the displacement of ionized and excited water molecules caused by charge migration. Oscillating charge migration is expected to happen a few femtoseconds after an energy transfer event, hence challenging its experimental observation given the short period of time [START_REF] Kuleff | Core Ionization Initiates Subfemtosecond Charge Migration in the Valence Shell of Molecules[END_REF]. Thus, RMS values have been empirically adjusted in other Monte Carlo track-structure codes due to the lack of corresponding measured data [START_REF] Kreipl | Time-and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation[END_REF][START_REF] Cobut | Monte Carlo simulation of fast electron and proton tracks in liquid water -I. Physical and physicochemical aspects[END_REF][START_REF] Tomita | Monte Carlo simulation of physicochemical processes of liquid water radiolysis[END_REF][START_REF] Uehara | Monte Carlo simulation of water radiolysis for low-energy charged particles[END_REF]. The RMS obtained in this work (1.24 nm), shorter than Geant4-DNA's default value (2 nm), increased the chance of contact reactions at the earliest times producing lower yields of • OH radical. For • OH radicals at 7 ps, this modification represented a reduction from 5.0 molecules per 100 eV (Ramos-Mendez et al., 2020) to 4.7 molecules per 100 eV (the measure data show 4.7±0.2 molecules per 100 eV [START_REF] Wang | Timedependent yield of the hydrated electron and the hydroxyl radical in D2O: A picosecond pulse radiolysis study[END_REF]). Geant4-DNA being a radiation transport code is not currently capable of simulating molecular dynamics of water molecules. Thus, effects like the binding of water molecules to biomolecules (DNA, proteins or lipids) that in principle affect the RMS displacement (see e.g., Shweta and Sen 2018), cannot be studied in detail at this point with our code. For e - aq , this parameter had a negligible effect as the products of auto-ionization that followed a dissociation event were handled by the one-step thermalization method [START_REF] Shin | Evaluation of the influence of physical and chemical parameters on water radiolysis simulations under MeV electron irradiation using Geant4-DNA[END_REF].

While it is a common practice, the comparison between time-dependent product yields calculated in pure liquid water and measured data under the presence of scavengers deviated from each other and should be avoided. In this work, two sets of data were selected to demonstrate this point. The temporal evolution of H 2 O 2 and H • reproduced the measured data within experimental errors only when the yields were calculated in the presence of scavengers. However, a more comprehensive evaluation of TOPAS-nBio under a wide range of scavengers, scavenger concentrations, and radiation qualities was outside the scope of this work and it is the subject of future work. A selection of radiobiologically relevant scavengers should precede that task given the scope of TOPAS-nBio.

On the other hand, for DNA damage simulations an encouraging agreement was found between calculated yields of SSB and DSB with measured data for low-LET radiation on a basepair level. A limitation of the applicability of IRT for smaller scales (e.g., including adenine, thymine, cytosine, and guanine nucleobases or atomic components) could be expected. Bluett and Green demonstrated [START_REF] Bluett | Competitive diffusion-influenced reaction of a reactive particle with two static sinks[END_REF]) that at such proximity, multiple reactive centers in the sugar-base system might occur, so the reaction rate between a reaction pair is affected by the presence of static neighboring species. The IRT method assumes that reactions between pairs occur in isolation, thus, the IRT method may give inaccurate results. Nevertheless, the use of observed reaction rates in this work provided a base-pair level resolution, which resulted in a general agreement within experimental uncertainties.

The estimated SSB efficiencies from the Monte Carlo data included a highly detailed plasmid DNA geometrical model. The estimated efficiency of 24% for SSB induction for • OH was within previously reported values between 24% to 44% obtained with cylindrical non-homogeneous kinetic model [START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF][START_REF] Udovičić | Yields of Single-Strand Breaks in Double-Stranded Calf Thymus DNA Irradiated in Aqueous Solution in the Presence of Oxygen and Scavengers[END_REF][START_REF] Klimczak | Irradiation of plasmid and phage DNA in water-alcohol mixtures: Strand breaks and lethal damage as a function of scavenger concentration Int[END_REF]. In that model, however, a homogenous cylindrical model representing a straight DNA segment was considered, which differed from the more detailed supercoiling DNA representation facilitated by TOPAS-nBio. For H • , the efficiency of 0.5% provided the best agreement at higher scavenger concentrations. This value was close to the 0.81% calculated by [START_REF] Aydogan | Monte Carlo Simulations of Site-Specific Radical Attack to DNA Bases[END_REF]. The differences are subtle and could be attributed to the different Monte Carlo codes used used by these authors. In their model, the authors used a straight linear DNA segment of 38 bp, with atomic resolution and multiple reactions sites on a base pair. However, our calculated efficiency was 12% higher than that measured experimentally [START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF]. This difference may be due to the limited resolution of our model, justified by the use of the variable rate coefficient for • OH + DNA, which disregards the individual reactions with nucleobases and other components. The reaction of • OH with nucleobases is in general several times faster than the reaction of • OH with deoxyribose [START_REF] Buxton | Critical Review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radicals in Aqueous Solution[END_REF]. Thus, by including reactions with DNA bases, an increment in the number of • OH reactions is expected, which might lead to a lower SSB efficiency. For that, reaction rate constants obtained in a B-DNA structure instead of an aqueous solution of DNA are needed. From the computational point of view, an effort of our group to find ways to extend the IRT method to handle DNA bases was reported in [START_REF] Tran | Assessment of DNA damage with an adapted independent reaction time approach implemented in Geant4-DNA for the simulation of diffusion-controlled reactions between radio-induced reactive species and a chromatin fiber[END_REF]. On the other hand, our estimation of the direct effect assumed an energy threshold for accumulated energy depositions of at least 17.5 eV, which produced 28% fewer SSB compared to, e.g., the 5-37.5 eV linear ramp threshold shown elsewhere (Zhu et al 2020b). A more suitable selection of the energy threshold must consider the transport of low energy electrons below 15 eV (avoiding the one-step thermalization model), which can induce DNA damage, as shown experimentally by [START_REF] Alizadeh | Biomolecular damage induced by ionizing radiation: The direct and indirect effects of low-energy electrons on DNA[END_REF]. In that regard, an extended physics package which improves the detail in the prechemical stage for handling Auger electrons and electron capture will be soon made available in Geant4 (Shin et al., 2021 under review). Experimentally, a dependence with the scavenging capacity of the SSB induction efficiency was further observed (Önal et al 1988). Even though we used a single constant value as a first approach, the assumptions made in this work were sufficient to reproduce the behavior of scavenging dependence of DSB yields measured for different plasmids.

Conclusions.

In this work, TOPAS-nBio was validated for simulating water radiolysis in liquid water and DNA damage for plasmids at low LET irradiation. Satisfactory agreement within experimental uncertainties was obtained, reconciling Monte Carlo calculations of water radiolysis from • OH and e - aq yields at the picosecond level. For pUC18 plasmids irradiated by 137 Cs γ-rays, calibration via the SSB induction efficiency was consistent with published efficiency values, leading to agreement with measure data within experimental uncertainties. TOPAS-nBio facilitated the implementation of experimental conditions of DNA irradiations including DNA geometry, source quality, prescribed dose, chemical parameters, and scavengers, exploiting the full potential of Geant4-DNA. As a result, an accurate, fast, and user-friendly Monte Carlo framework is provided by TOPAS-nBio/Geant4-DNA to evaluate DNA damage from first principles.

Figure 1

 1 Figure 1 Setup showing a two-stage simulation. The condensed-history MC simulation setup used to retrieve the secondary electron spectrum is shown on the left side. The track-structure MC simulation setup used to calculate SSB and DSB yields using supercoiled plasmid DNA is shown on the right side. Red lines correspond to few electron tracks. For more details, see the text.

Figure 2

 2 Figure 2 Time-dependent G values for fast electrons (1 molec./100 eV = 1.036x10 -7 mol J -1). TOPAS-nBio/Geant4-DNA simulated data: (solid line) pure liquid water calculations; (blue squares connected with dashed lines) simulations of scavenger systems for H 2 O 2 and H • as shown in table 1. Error bars represent statistical uncertainties, one standard deviation. Measured data: black and grey solid lines[START_REF] Ma | Scavenging the Water Cation in Concentrated Acidic Solutions[END_REF]; ☐[START_REF] Wang | Timedependent yield of the hydrated electron and the hydroxyl radical in D2O: A picosecond pulse radiolysis study[END_REF]; △[START_REF] Laverne | OH Radicals and Oxidizing Products in the Gamma Radiolysis of Water Source[END_REF]; ◼[START_REF] Bartels | Spur decay of the solvated electron in picosecond radiolysis measured with time-correlated absorption spectroscopy[END_REF]; ◇[START_REF] Shiraishi | Pulse-radiolysis study on the yield of hydrated electron at elevated temperatures[END_REF]; ✕ and ◆ (Pastina et al 1999); ○ (Hiroki et al 2002); ▲ (Huerta Parajon et al 2008).

Figure 3

 3 Figure 3 Total calculated SSB yields (direct plus indirect) as a function of • OH scavenging capacity (solid lines). Measured data for pUC18 is from ◼[START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF], ◆[START_REF] Milligan | Yield of single-strand breaks due to attack on DNA by Scavengerderived radicals[END_REF], and [START_REF] Milligan | Characterization of the reaction rate coefficient of DNA with the hydroxyl radical[END_REF]. Measured data for pUC19 ▲ is from[START_REF] Perry C C, Ramos-Méndez | Boronated Condensed DNA as a Heterochromatic Radiation Target Model Biomacromolecules acs[END_REF]. Percentage differences (☐) between calculated to measured data are shown in the bottom figure. The dotted lines are margins from experimental uncertainty. Error bars are displayed when bigger than the symbol, represent statistical uncertainties from Monte Carlo simulations, one standard deviation.

Figure 4

 4 Figure 4 Calculated G(SSB) as a function of the DNA concentration (◼ connected with solid line). Error bars, smaller than the symbols, represent statistical uncertainties, one standard deviation. Measured data is from[START_REF] Milligan | Variation of Single-Strand Break Yield with Scavenger Concentration for Plasmid DNA Irradiated in Aqueous Solution[END_REF]: pUC18 (○); pEC (△); and SV40 (◆).

Figure 5

 5 Figure 5 Calculated single and double-strand break yields as a function of hydroxyl radical scavenging capacity. Error bars represent statistical uncertainties, one standard deviation. Measured data: ○ pBR322[START_REF] Klimczak | Irradiation of plasmid and phage DNA in water-alcohol mixtures: Strand breaks and lethal damage as a function of scavenger concentration Int[END_REF], pBR322[START_REF] Tomita | Strand Break Formation in Plasmid DNA Irradiated in Aqueous Solution: Effect of Medium Temperature and Hydroxyl Radical Scavenger Concentration[END_REF].

Reactions for simulation of radiolysis in pure liquid water. Reactions with scavengers for DNA damage simulation.

  and (Huerta Parajon et al 2008)

			CH 3 OH + • OH  • CH 2 OH + H 2 O	9.7 x 10 8
			NO 3 -+ e aq - • NO 3 2-	9.7 x 10 9
			NO 3 -+ • H  • NO 3 2-	1.4 x 10 6
	Reactions with DNA	Reactions for calculating G(H • )
	Reaction	k obs (/M/s)	Reaction	k obs (/M/s)
	• OH + DNA	Variable	•H + HCO 2 - H 2 +•CO 2 -	2.1 x 10 8
	H • + DNA	0.03 x 10 9	e -aq + NO 3 - NO 3 2-	9.7 x 10 9
	e -aq + DNA	0.01 x 10 9	Br + • OH 	1.1 x 10 10
	Reaction	k obs (/M/s)
	e -aq + e -aq  H 2 + OH -	5.5 x 10 9	Reaction	k obs (/M/s)
	e -aq + H 3 O +  H •	2.3 x 10 10	e -aq + O 2  O 2 -	1.9 x 10 10
	e -aq + H •  H 2 + OH -	2.5 x 10 10	H • + O 2  HO 2	2.1 x 10 10
	e -aq + • OH  OH -	3.0 x 10 10	• OH + DMSO	7.1 x 10 9
	e -aq + H 2 O 2  OH -+ • OH	1.1 x 10 10	H • + DMSO	2.7 x 10 6
	H 3 O + + OH - H 2 O	14.3 x 10 10	e -aq + DMSO	3.8 x 10 6
	H • + H •  H 2	7.8 x 10 9
	H • + • OH  H 2 O	1.55 x 10 10
	H • + H 2 O 2  • OH + H 2 O	9.0 x 10 7	Reactions for calculating G(H 2 O 2 )
	• OH + • OH  H 2 O 2	5.5 x 10 9	Reaction	k obs (/M/s)

Table 2

 2 Displacements of physical and dissociation products implemented in Geant4-DNA.

	Physical product Dissociation	𝑅𝑀𝑆 𝐻 2 𝑂 +	or	Displacement
		products	𝑅𝑀𝑆 𝐻 2 𝑂 *
	H 2 O +	H 3 O + + • OH (1)	2.0 nm		H 3 O +	0 or 0.8 nm (2)
					• OH	0.8 or 0 nm (2)
	H 2 O *	H • + • OH	2.4 nm		H •	17/18 𝑅𝑀𝑆 𝐻 2 𝑂 *
					• OH	1/18 𝑅𝑀𝑆 𝐻 2 𝑂 *
		H 2 + 2 • OH	0.8 nm		H 2	2/18 𝑅𝑀𝑆 𝐻 2 𝑂 *
					• OH • OH	16/18 16/18	𝑅𝑀𝑆 𝐻 2 𝑂 * 𝑅𝑀𝑆 𝐻 2 𝑂 *	+ 0.55 nm + 0.55 nm (3)
		H 2 + • OH + OH -0.8 nm		H 2	2/18 𝑅𝑀𝑆 𝐻 2 𝑂 *

Acknowledgments

This work was supported by NIH/NCI R01 CA187003 (TOPAS-nBio). V.Š. and K.S. were supported by the Czech Science Foundation project no. 17-03403Y. N.D-K. is a doctoral student from Programa de Doctorado en Ciencias Física Aplicada, Benemérita Universidad Autónoma de Puebla and received fellowship 2019-000002-01NACF-05024 from CONACYT, México.

Faddegon B, Ramos-Méndez J, Schuemann J, McNamara A, Shin J, Perl J and Paganetti H 2020 The TOPAS tool for particle simulation, a Monte Carlo simulation tool for physics, biology and clinical research Phys. Medica 72 114-21 Online: https://doi.org/10.1016/j.ejmp.2020.03.019 Friedland W, Schmitt E, Kundrát P, Dingfelder M, Baiocco G, Barbieri S and Ottolenghi A 2017 Comprehensive track-structure based evaluation of DNA damage by light ions from radiotherapy-relevant energies down to stopping Sci. Rep. 7 45161 Online: http://www.nature.com/articles/srep45161 Fulford J, Nikjoo H, Goodhead D T and O'Neill P 2001 Yields of SSB and DSB induced in DNA by Al K ultrasoft X-rays and α-particles: comparison of experimental and simulated yields Int. J. Radiat. Biol. 77 1053-66 Online: http://informahealthcare.com/doi/abs/10.1080/09553000110069308 Green N J B, Pilling M J and Clifford P 1990 Stochastic Modeling of Fast Kinetics in a Radiation Track Society 94 251-8 Online: http://pubs.acs.org/doi/abs/10.1021/j100364a041 Hiroki A, Pimblott S M and Laverne J A 2002 Hydrogen peroxide production in the radiolysis of water with high radical scavenger concentrations J. Phys. Chem. A 106 9352-8 Huang J, Schlick T and Vologodskii A 2001 Dynamics of site juxtaposition in supercoiled DNA Proc. Natl. Acad. Sci. U. S. A. 98 968-73 Huerta Parajón M 2010 HYDROGEN ATOM FORMATION IN THE GAMMA AND HEAVY ION RADIOLYSIS OF AQUEOUS SYSTEMS A thesis submitted to The University of Manchester for the (The University of Manchester) Huerta Parajon M, Rajesh P, Mu T, Pimblott S M and LaVerne J A 2008 H atom yields in the radiolysis of water Radiat. Phys. Chem. 77 1203-7 Incerti S, Baldacchino G, Bernal M, Capra R, Champion C, Francis Z, GuÈye P, Mantero A, Mascialino B, Moretto P, Nieminen P, Villagrasa C and Zacharatou C 2010a THE Geant4-