Bringing in vitro analysis closer to in vivo: Studying doxorubicin toxicity and associated mechanisms in 3D human microtissues with PBPK-based dose modelling.

Doxorubicin (DOX) is a chemotherapeutic agent of which the medical use is limited due to cardiotoxicity. While acute cardiotoxicity is reversible, chronic cardiotoxicity is persistent or progressive, dose-dependent and irreversible. While DOX mechanisms of action are not fully understood yet, 3 toxicity processes are known to occur in vivo: cardiomyocyte dysfunction, mitochondrial dysfunction and cell death. We present an in vitro experimental design aimed at detecting DOX-induced cardiotoxicity by obtaining a global view of the induced molecular mechanisms through RNA-sequencing. To better reflect the in vivo situation, human 3D cardiac microtissues were exposed to physiologically-based pharmacokinetic (PBPK) relevant doses of DOX for 2 weeks. We analysed a therapeutic and a toxic dosing profile. Transcriptomics analysis revealed significant gene expression changes in pathways related to "striated muscle contraction" and "respiratory electron transport", thus suggesting mitochondrial dysfunction as an underlying mechanism for cardiotoxicity. Furthermore, expression changes in mitochondrial processes differed significantly between the doses. Therapeutic dose reflects processes resembling the phenotype of delayed chronic cardiotoxicity, while toxic doses resembled acute cardiotoxicity. Overall, these results demonstrate the capability of our innovative in vitro approach to detect the three known mechanisms of DOX leading to toxicity, thus suggesting its potential relevance for reflecting the patient situation. Our study also demonstrated the importance of applying physiologically relevant doses during toxicological research, since mechanisms of acute and chronic toxicity differ.

A C C E P T E D M A N U S C R I P T processes resembling the phenotype of delayed chronic cardiotoxicity, while toxic doses resembled acute cardiotoxicity. Overall, these results demonstrate the capability of our innovative in vitro approach to detect the three known mechanisms of DOX leading to toxicity, thus suggesting its potential relevance for reflecting the patient situation. Our study also demonstrated the importance of applying physiologically relevant doses during toxicological research, since mechanisms of acute and chronic toxicity differ.
A C C E P T E D M A N U S C R I P T

Introduction
Doxorubicin (DOX, also known by the brand name Adriamycin) is a chemotherapeutic agent belonging to the class of anthracyclines, which are cytostatic antibiotics. DOX was isolated from Streptomyces peucetius var. caesius in 1967 [1] and it is still one of the most frequently used anti-cancer agents for treating a variety of cancers (e.g. hematological cancers, soft-tissue tumors, and solid tumors), even though it causes severe side effects. Next to nausea, vomiting, alopecia, myelosuppression, stomatitis, and gastrointestinal disturbances, DOX is known to induce cardiotoxicity. Cardiotoxicity is classified as acute or chronic. Acute cardiotoxicity may already occur during treatment with a single high dose or within 2-3 days after multiple DOX treatments [2,3]. Though its prevalence is higher than for chronic cardiotoxicity, these side effects are reversible and clinically manageable. A chronic cardiotoxic phenotype may emerge between one month and decades after DOX treatment. This persistent or progressive cardiotoxicity is dose-dependent and irreversible. In cases where congestive heart failure has developed, the prognosis for the patient is poor, with a mortality rate of 50% within one year [3][4][5][6][7].
Given this dose-dependency, a maximum cumulative dose of 450 -600 mg/m 2 is recommended for DOX [8,9]. DOX mechanisms of action are not fully understood yet. Nevertheless, it is generally accepted that the main mode of action is related to killing dividing cells. These effects are more severe for tumor cells since these cells divide more rapidly than non-cancer cells. However, DOX is not tumor cell-specific and can accumulate in the nucleus and mitochondria of heart, liver and blood cells, thereby contributing to toxic side effects, in particular chronic cardiotoxicity, even at therapeutic doses [5]. DOX has been found to 1) intercalate into DNA, 2) target DNA topoisomerases, and 3) generate reactive oxygen species (ROS) [10][11][12]. The first two processes inhibit unwinding of DNA, DNA replication, RNA transcription and protein biosynthesis. As a result, proliferation of dividing cells is also inhibited [9,13]. This is thought to be the efficacy of the anti-cancer effects of DOX, while ROS generation is mainly ascribed to toxic effects [14]. ROS can be generated as a result of mitochondrial dysfunction, but also by the oxidative semiquinone formed at complex I of the electron transport chain (ETC) during DOX metabolism. The induced oxidative stress may damage cells and cause cell death [8,11]. Other detrimental actions of DOX can be related to death of non-cancer cells or to decreased cardiomyocyte functioning, which may partly result from mitochondrial dysfunction causing an imbalance in cellular energetics. Therefore, any effect of DOX on mitochondria structure or function can cause cardiomyocyte dysfunction and thereby cardiotoxicity [15]. The detrimental actions of DOX, known to occur in vivo, can thus be summarized into 3 toxicity processes: cardiomyocyte dysfunction, mitochondrial dysfunction and cell death. However, it should be noted that the split between efficacy and detrimental actions is not fully justified because of overlapping processes, such as oxidative stress and cell death.
To predict molecular mechanisms underlying long term toxicity, toxicological risk assessment traditionally relies on animal models [16]. However, notably due to increasing ethical pressure, the field has to reduce the amount of animal experiments. In this article, we present an innovative in vitro experimental design aimed to detect cardiotoxicity by obtaining a complete view of the induced mechanisms of a compound. By better reflecting the human in vivo conditions, we aim to find an in vitro model capable of reliably replacing animal models used within the field of toxicology and drug safety testing. First, instead of the regularly used monolayer cell cultures, a Human 3D cardiac microtissue (InSphero, SWL) model was used, which better resembles the in vivo heart [17][18][19]. This in vitro spheroid cell model contained approximately 4000 iPSC-derived human cardiomyocytes and 1000 cardiac fibroblasts. Second, in contrast to using a relatively high dose to treat the cells, physiologically based pharmacokinetic (PBPK) modelling [20] was included within the study design to better reflect the in vivo exposure. PBPK simulates the absorption, distribution, metabolism and excretion of a specific dose of DOX within the human body, therefore enabling us to predict the concentration to which a specific organ is exposed over time. The PBPK model was used to design a two weeks repetitive dosing profile in which microtissues were exposed to a decreasing dose over each day by means of 3 medium changes.
The microtissues were exposed to either a therapeutic or a toxic dose of DOX using this PBPK-based repetitive dosing profile, and functional parameters for mitochondrial function and programmed cell death were assessed. Finally, by applying next-generation total RNA-sequencing to ribo-depleted RNA samples, we were able to analyze global changes in gene expression. By using this innovative experimental in vitro experimental design, we were able to successfully detect biological changes in the three main identified toxicity processes of DOX (cardiomyocyte dysfunction, mitochondrial dysfunction and cell death).

PBPK model establishment
A PBPK model for DOX was established with the PBPK modelling software PK-Sim following a previously described workflow [21]. In particular, the model was used to quantify interstitial heart concentrations following administration of therapeutic and toxic drug doses. These concentration profiles were subsequently translated into discrete daily exposure profiles (Figure 1a, b). During model establishment, model predictions of blood plasma DOX concentrations were compared with literature data for parameter identification [22]. In a next step, the model was validated with further PK data4 ( Figure 1c). Simulations of DOX interstitial heart concentration profiles following a therapeutic (green) and toxic (red) once daily dosing regimen with experimental sampling points (black circles). b) Exemplary translation of PBPK simulations into discrete in vitro assay concentration with equal exposure, applied to cardiac microtissues for the first three days. c) Validation of PBPK model predicted DOX blood plasma concentrations with literature data [23].

Exposure
To accomplish the PBPK-based repetitive dosing profile described above, the medium of the microtissues was changed three times daily on working days with the PBPK-determined DOX concentrations (Sigma D1515, Cat #25316-40-9). These were administered using stock solutions of the compound dissolved in DMSO. As control, microtissues exposed to similar end concentrations of DMSO were used. 7 time points (2,8,24, 72, 168, 240 and 336 hours) were sampled in triplicate during the 2 weeks treatment period. For each sample, 36 microtissues were incubated separately, and subsequently pooled before RNA extraction in order to generate sufficient input material for NGS analysis.

RNA sequencing & Data processing
From the exposed microtissues, total RNA was isolated using Qiagen AllPrep DNA/RNA/miRNA Universal Kit (Cat #80224). The sample was depleted of ribosomal RNA using the Illumina RiboZero Gold kit (Cat #MRZG12324) and prepared for sequencing using Lexogen SENSE total RNA library preparation kit (Cat #009.96). After library preparation, the samples were sequenced on the HiSeq2500 (100bp paired-end).
A pool of all DOX samples was sequenced on all 8 lanes of a flow cell. The same was done for the control samples. Control samples were sequenced with an average of 51.2 million raw reads, therapeutic DOX with 20.3 million reads and toxic DOX with 18.7 million. Because the sequencing reads still contained Lexogen adapter sequences, the first 12 bases of the 5'end of all reads were removed using Trimmomatic version 0.32 [24]. The quality of the sequencing data was checked using FastQC version 0.11.3 [25] before and after trimming.

Gene expression analysis
For each sample, sequencing reads were mapped to the GENCODE version 26 (GRCh38) -Ensembl 88 (Comprehensive gene annotation) using the Genedata Profiler TM version 10.1, which incorporates the splice junction mapper STAR version 2.5.2b and Cufflinks based algorithms for transcript abundance estimation.
The quality of the samples was assessed according to the amount of (mapped) reads, Cook's distance, hierarchical clustering, principal component analysis, and sample dispersion. 1 replicate of DOX toxic 240h was excluded because of low read count.
Differentially expressed genes (DEGs) were determined using DESeq2 version 1.14.1 [26]. For each specific time point, 2 comparisons were performed: 1) therapeutic dose vs DMSO control; 2) toxic dose vs DMSO control. Samples of all time points for both doses were normalized together and the contrast function was used to extract differentially expressed genes for each time point, ensuring comparability between time points. Genes were being considered as differentially expressed when q<0.05.

Biological interpretation
For each dose and time point, lists of DEGs (depleted of ribosomal genes) were used as input for pathway over-representation analysis (ORA) using CPDB release 32 [27]. The Reactome database version 61 was selected because it is curated, focuses on molecular mechanisms and includes pathway hierarchy [28]. Results from this analysis was visualized using Omix visualization version 1.9.20 [29].

Functional assessments: ATP measurement, Caspase3/7 assay and Oxygen consumption
For measurement of ATP, Promega's CellTiter Glo 3D (Cat #G9683) was used according to manufacturer's protocol. In short, the microtissues were incubated for 30min with luciferase reagent and the luminescence was measured.
Mitochondrial function after 2 hours and 7 days of DOX treatment was assessed by measurement of extracellular oxygen consumption using Luxcel's MitoXpress ® Xtra Oxygen Consumption Assay (Cat #MX- according to manufacturer's protocol. In short, oxygen quenches the MitoXpress ® Xtra probe, making the measured fluorescent signal, inversely proportional to the oxygen concentration.

Mitochondrial model
The computational mitochondrial model [30] for proteomics analysis was extended to evaluate the expression changes of genes functioning in mitochondria. The original model contained reactions of metabolites and involved protein complexes. For each protein complex represented in the model, the corresponding genes coding for each individual protein subunit were added.

Results
In order to investigate mechanisms of DOX-induced toxicity, human 3D cardiac microtissues were exposed to DOX for two weeks with either a therapeutic dose (DTHE) or a toxic dose (DTOX) using a PBPK based dosing profile. To this end, the medium of the microtissues was changed three times per workday with a decreasing dose. Samples were generated after 2, 8, 72, 168, 240 and 336 hours of DOX exposure for both doses (Figure 2A).
Molecular effects of DOX on the cardiac cell were investigated using RNA-sequencing. The amount of differentially expressed genes (DEGs) increased over time ( Figure 2B), except for a temporary decrease in DEGs at 240 hours of DTHE. Consistently across all time points, more DEGs were detected during exposure with DTOX than with DTHE. Of these DEGs, 80% were protein-coding genes.

Identifying known DOX pathways in vitro
An overview of the effects of DOX on the whole cardiac cell was investigated to identify processes linked to DOX toxicity. This was done through pathway overrepresentation analysis (ORA). Overrepresented pathways were visualized in supplementary figure 1. More significantly overrepresented pathways were identified for DTHE compared to DTOX. This in combination with the higher amount of DEGs detected for DTOX compared to DTHE indicates less specific cellular reactions for DTOX. While therapeutic DOX exposure affects many cellular processes, the most significantly affected pathways throughout DOX treatment are "TCA cycle & respiratory electron transport" (qmedian 2.92•10 -6 ) and "Striated muscle contraction" (qmedian 9.3•10 -6 ). These pathways can be related back to detrimental actions of DOX. Therefore, toxicity processes (cardiomyocyte dysfunction, mitochondrial dysfunction and cell death) were investigated in   Both DOX doses clearly induced "striated muscle contraction", the main pathway related to cardiomyocyte dysfunction. Both doses displayed also a decreasing amount of upregulated genes over time combined with an increase of downregulated genes. For DTOX, the amount of downregulated genes was always higher than of upregulated genes, while DTHE displayed more upregulated genes until 240 hour of treatment. In "cellular senescence", more DEGs were observed for DTOX compared to DTHE, while the pathway was found significantly overrepresented for DTHE at more time points than for DTOX. This phenomenon is common in ORA and caused by the higher amount of DEGs and less specific cellular reactions of DTOX. More striking were the differences between doses with respect to metabolic processes within the mitochondria. While DTHE induced upregulations over time, DTOX induced downregulations at the end of the DOX treatment period(168h -336h). Because mitochondria are involved in cell death and ATP generation necessary for cardiac contractile function, this dose-dependent difference in DOXinduced expressions of genes involved in mitochondrial processes may represent the underlying cause for differences in phenotypes for chronic or acute cardiotoxicity. Therefore, it is of interest to analyze mitochondrial dysfunction in more detail.

Functional assessment of microtissue viability
With respect to the ATP content of DOX-treated cells, the two analyzed doses ( Figure 4A

Mitochondrial effects of DOX exposure
To elucidate the effects of DOX on genes functioning in mitochondria, DEGs were mapped onto the computational mitochondrial model [30] adapted for transcriptomics data ( Figure 5). As expected, the mitochondrial model showed major differences between the two doses. DTHE displayed a majority of upregulated genes, while DTOX displayed more downregulations, especially at later time points.   Figure 5) was activated to the same extent for both doses, though there are slight differences in the final step of ROS detoxification, in which H2O2 is converted to water.

Detoxification of ROS Detoxification of ROS (bottom left of
While for DTHE this reaction was mainly catalyzed by peroxiredoxin (PRX), for DTOX this reaction was carried out by both peroxiredoxin (PRX) and glutathione peroxidase (GPx). Overall, the upregulations found in these ROS detoxification genes indirectly indicated an increase in ROS formation in response to DOX. Dose-dependent differences were most noticeable in the complexes of the ETC (top left of Figure 5). The overall amount of genes affected (minimally at 1 time point) was similar between the doses, with complex I expressing the most DEGs. Because of the difference in amount of genes per complex, comparison in percentage of complex gave a better overview of the DOX effect (Table 1). Complex III was mostly affected for DTHE, while complex V (also known as ATP synthase) was highly affected by both doses. Complex II was the least affected by DOX exposure. Dose-dependent differences were mostly observed through the direction of expression changes. Affected genes for DTHE were mostly upregulated, while DTOX displayed an increase in downregulated genes, especially in later time points. Names of genes that were most highly affected by DOX exposure (significantly affected minimally at 4 time points) were also included in Table 1. Figure 6 depicts the log2 fold changes of these genes. Only gene expression changes of highly affected complex V genes were highly correlated to the changes observed in ATP Overall, DOX-induced gene expression changes in the ETC were observed at both doses, though these changes indicate different underlying processes between the doses, discriminating between the phenotype of delayed chronic cardiac toxicity and acute cardiotoxicity.

Discussion
The main objective of this study was to validate our novel in vitro experimental design for its use in predictive toxicological research. For this purpose, we analyzed cardiotoxic effects of doxorubicin (DOX).
Our in vitro approach was designed to improve the resemblance to in vivo situations. This was achieved by using novel human 3D cardiac microtissues. This spheroid cell model contains a co-culture of iPSCderived human cardiomyocytes and cardiac fibroblasts, better resembling the in vivo heart than monoculture monolayers. Furthermore, while standard toxicological research is being done with short time series and exposures at one specific dose, we applied a 2 week PBPK-based repetitive dosing profile that resembles in vivo drug concentrations in patients with one drug administration per day. The dose of DOX to which a patient's heart is exposed was calculated based on pharmacokinetic simulations of the absorption, distribution, metabolism and excretion of DOX. In vitro, the medium of the cardiac microtissues was changed three times daily with calculated DOX doses, resulting in the highest dose in the morning and the lowest dose in the evening of each working day during the 2 week exposure. The experimental design also included the exposure of two different dosing profiles, one based on therapeutic doses, as administered to a patient, and one based on a toxic dose. From these samples, RNA was extracted and sequenced to determine if this novel in vitro approach is able to retrieve the in vivo toxicity mechanisms known for DOX.
We used three detrimental actions of DOX which are known to occur in vivo (cardiomyocyte dysfunction, mitochondrial dysfunction and cell death) to determine that the in vitro approach was able to retrieve the in vivo toxicity pathways known for DOX [8,11,14,15]. Results from ORA indicated less specific cellular reactions for DTOX, observed by the combination of more DEGs and less overrepresented pathways compared to DTHE. Furthermore, ORA identified toxicity of DOX by gene expression changes in the pathways of "striated muscle contraction", "TCA cycle & respiratory electron transport" and "cellular senescence". Gene expression changes in "striated muscle contraction" (cardiac function) and "cellular senescence" (cell death) displayed similar patterns, while mitochondrial processes ("TCA cycle & respiratory electron transport") displayed significantly upregulated processes for DTHE and non- complexes. This process is also supported by the expression profile in mitochondrial encoded genes for subunits of complex I (MT-ND genes), which displayed changes in toxic dose, but not in therapeutic dose.
The downregulations in gene expression at the later time points of toxic DOX exposure indicates the failure of ATP production through ETC. The decrease in ATP content of the microtissues was a continuous process throughout the exposure, a pattern that was also observed for ATP5S, the most affected gene in complex V. This gene, also known as factor B, is responsible for the energy transduction activity of the ATP synthase complex [33]. Though, it is unlikely that a defect of this one gene was solely responsible for the ATP decrease. Decreased effectiveness of the ETC may also be due to decreased oxygen availability, oxygen being consumed by DOX during ROS generation. Therefore, the ETC has less protons available to establish the proton gradient necessary for ATP production [14]. However, it seems more likely that the heterogeneity of DOX exposure causes malfunctioning cells in the outer layers of the microtissue to die, leaving a decreased amount of viable cells at the core to produce ATP. Because In microtissues exposed to therapeutic doses, the DOX effect was also not equally distributed over the mitochondria present in the cells of the cardiac microtissue. Only a small amount of cells in the outer layer of the microtissue died, as could be seen from the small increase in caspase activity and a decrease in microtissue size that is hardly visible. However, an increase in gene expression of mitochondrial and ETC genes was necessary in order to maintain cellular ATP content, indicating a possible adaptation response in which less affected mitochondria compensate for highly affected mitochondria. However, this adaptation response may not be sufficient for maintaining cardiac function, because genes involved in striated muscle contraction were already showing decreasing expression level. Since the regenerative capacity of the heart is limited to ventricular remodelling, cumulative toxicity surpasses the point to which the organ can adapt, thereafter resulting in a delayed phenotype of heart failure, as is observed in patients with chronic cardiotoxicity after DOX treatment [4,15].

Conclusion
Overall, these results confirm that our in vitro approach, based on 3D human microtissues exposed to an innovative 2 weeks PBPK dosing profile of DOX, is able to retrieve the known mechanisms of DOX toxicity. This new in vitro toxicity model can be extended to other 3D human microtissues, such as the already available liver model from primary hepatocytes.
The insights from the current study demonstrate the importance of applying physiologically relevant dosing profiles during toxicological research. Many studies on DOX effects have been performed with single and/or extremely high doses that do not reflect -and may even obscure -the true mechanisms of toxicity. Therefore, controversy in published results must be treated cautiously because of the risk of oversimplification and mixture of the different mechanisms of acute and chronic effects [14,34].
Ultimately, we believe that this novel approach, better mimicking the human in vivo condition, could replace animal models within the field of toxicology and drug safety testing.

Conflict of Interest
The authors declare that there are no conflicts of interest.