A Dynamic Model for Estimating the Interaction of ROS–PUFA–Antioxidants in Rabbit

Defining optimal nutrition in animals and humans remains a main scientific challenge. The objective of the work was to develop a dynamic model of reactive oxygen species (ROS)–polyunsaturated fatty acid (PUFA)–antioxidant homeostasis using the rabbit as a model. The problem entity was to evaluate the main metabolites generated from interactions between traits included in the conceptual model and identified by three main sub–models: (i) ROS generation, (ii) PUFA oxidation and (iii) antioxidant defence. A mathematical model (VENSIM software) that consisted of molecular stocks (INPUTs, OUTPUTs), exchange flows (intermediate OUTPUTs) and process rates was developed. The calibration was performed by using standard experimental data (Experiment 1), whereas the validation was carried out in Experiments 2 and 3 by using supra–nutritional dietary inputs (VIT E+ and PUFA+). The accuracy of the models was measured using 95% confidence intervals. Analytical OUTPUTs (ROS, PUFA, Vit E, Ascorbic acid, Iso–/NeuroProstanes, Aldehydes) were well described by the standard model. There was also good accuracy for the VIT E+ scenario, whereas some compensatory rates (Kc1–Kc4) were added to assess body compensation when high levels of dietary PUFA were administered (Experiment 3). In conclusion, the model can be very useful for predicting the effects of dietary treatments on the redox homeostasis of rabbits.


Introduction
Reactive oxygen species (ROS) such as peroxides, the superoxide anion, the hydroxyl radical and singlet oxygen have complex interactions with the different classes of biomolecules and physiological processes of living organisms [1], thus playing fundamental roles in health and disease [2].
In isolated mitochondria, it is estimated that 0.1-2% of oxygen consumed produces ROS [3]. Depending on their relative concentrations and site of production, these molecular species can damage different cellular targets, including proteins [4], nucleic acids [5] the primary defences are not sufficient, other mechanisms, such as transcription of genes and transduction of essential proteins, interfere with the oxidative response (i.e., the cytochrome P450 [CYP450] family [32]). Furthermore, the tocopheroxyl radical can be reduced back to tocopherol by the activity of a redox antioxidant system that includes AA and thiol-dependent enzymes (e.g., glutathione reductase) as the main players [12], and the scavenging function of vitamin E (Vit E) cooperates with that of the membrane peroxidase, especially phospholipid hydroperoxide 4 (GPx4) [33].
Quantifying ROS, and, in general, the oxidation products of lipids and other biomolecules, remain very challenging. Current analytical tools provide different opportunities, but most of these exhibit a lack of specificity and sensitivity. Most common measurements are based on oxidation biomarkers, including aldehydes and isoprostanoids [34][35][36]. However, using these biomarkers without evaluating the molecules involved in the oxidative processes of tissues provides little information. Moreover, they do not represent the entire metabolic system, but only a portion of the mechanism.
The objective of this work was to build a model to simulate redox homeostasis of the rabbit in which the dynamic interactions of blood plasma indicators of ROS production and variations in PUFA and antioxidant levels (i.e., Vit E and AA) are considered during the response to metabolic and environmental stimuli (such as dietary intake of Vit E and PUFAs).

Cognitive Map Used for Definition, Calibration and Validation of the Model
The method reported by Sargent [37] was applied to build the model. Accordingly, we defined the Problem Entity, as well as the Conceptual and the Structural (i.e., computerised) models (Sections 2.1.1-2.1.3), integrating it with the states of Trucano et al. [38] and Dimokas et al. [39] to define the verification and validation procedures (Section 2.2).
Once we verified that the Structural model represents the Conceptual Description, we calibrated the model (showing the degree to which the model accurately fit the in vivo values (experiment 1, see Section 2.3) and successively validated it with another set of data (experiments 2 and 3, see Section 2.3).

Problem Entity
As introduced earlier, evaluation of single products generated during the redox homeostasis of a system (i.e., aldehydes, isoprostanoids and antioxidants) from the interaction of multiple factors (ROS-PUFA-antioxidants) cannot describe the complexity and dynamic interactions of internal and external factors. Figure 1 shows the map of the conceptual model used in this study. We considered three main sub-models for simulations and validation of the main biological interactions in the redox-homeostasis system of the rabbit plasma: (i) ROS generation, (ii) PUFA oxidation and (iii) antioxidant defences. The rate of every reaction depends on the metabolic balance of the model (see Section 2.1.3).

Structural Model
Inflows and outflows represent 'gaining' and 'losing' molecules, respectively. The 'dynamic' flow of substrates and products of a reaction, including transient intermediates (i.e., hydroperoxides, peroxides, Vit E oxidation, etc.) is represented by the reaction rates as reported in the following sections. The time considered is 1 h (3600 s), which represents about one fifth of the time for standard feed digestion in the rabbit (4-6 h) [36]. To build the operation map in VENSIM software [40], different traits were analytically determined (Section 2.5) and classified as INPUTs and OUTPUTs ( Figure 2 [12]. The structure of the dynamic model built with the VENSIM software. The model is split into three main sub-models: ROS generation (blue), PUFA oxidation (orange) and antioxidant defence (yellow) (as detailed in Figure 1). The variables inside the angle brackets (<>) are reported as a shadow of the stock variables. Abbreviations: AA, ascorbic acid; AA ox, oxidised ascorbic acid; PUFA, polyunsaturated fatty acids; ROS, reactive oxygen species; Vit E, vitamin E; Vit E ox, oxidised vitamin E.

Sub-Model 1: ROS Production
For ROS, the physiological concentrations and basal rate of generation are described in the model system as ROS input and kROS, respectively. These variables are assumed to represent both the standard metabolism of a living system (i.e., the respiratory process, movement and cell metabolism [16]) and the compensatory response to external stimuli (Kc4 and Kc4.1; e.g., hyperventilation, nutrition, stress, environmental perturbation, etc.), which can increase and even exceed the antioxidant capacity of the system and lead to oxidative stress [41]. ROS output (Equation (1)) is the sum of ROS input and ROS generated through the neutralisation of Aldehydes (A6; Equation (11)). Furthermore, a numeric value (i.e., 400) has been added to convert the units of the considered variables to moles. A1 (Equation (2)) represents the reactive molecules originated by the reaction steps of Vit E.
Sub-Model 2: Oxidative Thrust In the model, PUFA represents the main substrate of oxidation. PUFA can be introduced with the diet (input PUFA) to follow different metabolic routes associated with specific metabolic pathways (A1-A6; Equations (3)-(12)).
The third pathway is β-oxidation (βOxidation; Equation (5)) [42] that, independently of ROS activity, promotes LC-PUFA catabolism (reported in the system as kUNS) to produce energy in the mitochondria with a rate (kβOx) affected by the energy balance of the different tissues, negatively correlating with storage of PUFA in cells [43].
The fourth route is autoxidation of PUFA to produce isoprostanoids (IsoNeuroPs). These are generated from LC-PUFA in cell membranes. After oxidation, they are cleaved by phospholipase and released in the interstitial fluids. Consequently, IsoNeuroPs are found in the plasma and can be eliminated through phospholipase A2 [44] with urine (IsoNeuro excretion; Equation (6)) at a specific excretion rate (kIsoNeuro excretion). Their concentrations and autoxidation (AutoOx; Equation (7)) are assumed to depend on PUFA precursors and/or ROS (kIsoNeuroPs). Compensations introduced for this pathway are Kc1 and Kc1.1 (see Section 2.2).
The fifth pathway is PUFA oxidation that correlates with the unsaturation level (kUNS) of FA, ROS and antioxidant activity [45]. PUFA oxidation produces Peroxides (initialisation step, A3; Equation (8)), a phenomenon that can be prevented, at least in part, by the H atomdonating activity of Vit E. This is also expected to influence Hydroperoxide levels, depending on the downstream activity of membrane peroxidases (primary metabolites resulting from termination activity, A4; Equation (9)). If not properly scavenged, lipoperoxyl radicals and lipoperoxides undergo extensive fragmentation to form meta-stable-molecules, such as Aldehydes (secondary metabolites, A6; Equation (11)). In particular, one molecule of hydroperoxide is broken down into two molecules of aldehydes (final Aldehydes) [46]. The Aldehyde generation depends on antioxidants (reported as kAldehydes), and on molecules' susceptibility to oxidation (kPerox, A5; Equation (10)). Finally, aldehydes can be excreted by urine (Aldehyde excretion), assuming that the rate of excretion is related to their concentrations (kAldehyde excretion; Equation (12)). This rate can change in relation to the activation of compensatory activities aimed at reducing the damage to lipids and proteins (Kc2 and Kc2.1; as reported in the following scenarios Section 2.3). It should be underlined that the concentrations of these latter molecules could be affected by mechanisms of re-absorption, transport and/or excretion that are difficult to estimate with our experimental plan.
Note that 1000 is the scale factor to compare Vit E to peroxides.

Sub-Model 3: Antioxidant Defence
Antioxidants are substances that delay or prevent the oxidation of other molecules. In our model, Vit E is the main fat-soluble antioxidant capable of preventing PUFA peroxidation. It works both as a scavenger of peroxyl radicals and an inhibitor of oxidoreductases involved in the enzymatic oxidation of PUFA. The antioxidant activity of Vit E results in the formation of α-TQ (Vit Eox and Vit Eox remained or Vit Eox transit, A7; Equations (13) and (14), respectively) with a molar rate of 1/2 with respect to peroxyl radicals scavenged during the H atom transfer process. Vit E in the model changes according to different metabolic pathways (A7-A9) and together with AA is regulated by Equations (13)- (20).
The availability of Vit E is modulated by a series of exogenous and endogenous factors [47]. Under physiological conditions, the main exogenous factor is dietary intake (Vit E input) [12]. Furthermore, Vit E (A8; Equation (15)) stored in the cells or other depots (cell storage) and recovered by recycling tocopheryl radicals formed during lipoperoxyl radical scavenging (Vit E recycling; Equation (16)) is the endogenous counterpart acting on the inner pool of the vitamin. Indeed, Vit E is regenerated by the intervention of AA, which acts as a secondary antioxidant (A9; Equation (17)), with a specific rate indicated as kVit E_AA.
Vit E not recycled after oxidation or in excess is assumed to be excreted in the urine or bile (Vit Eox excretion) with a concentration-dependent rate (kVit Eox excretion; Equation (18)). Compensatory mechanisms considered in this respect are indicated as Kc3 and Kc3.1, which include CYP450 detoxification and other excretion processes [48]).
The same recycling process is considered in the case of AA that is oxidised (AAox) according to AAox transit. Thus, it is restored back to the reduced form (AA recycling; Equation (19)) by trans-hydrogenases modulated by glutathione reductase (GR) [49] and specific ascorbate reductases [50], with a rate identified in the model as kAA enzyme. Similarly to Vit E, AA is modulated by dietary intake (AA input), AAox excretion (Equation (20)), the rate of excretion (kAAox excretion) and AA recycle.

Verification and Validation of the Model
We calibrated the model by using the data from Experiment 1 (standard model) to evaluate its accuracy. We validated the model by using data from Experiments 2 and 3, where the main INPUTs were changed (VIT E+ and PUFA+ scenarios; Figure 4); the fitness of the results was assessed.
We confirmed the model if the model outcomes of Experiments 2 and 3 were within the 95% CI (e.g., VIT E+). Otherwise, starting from the examination of final and intermediate OUTPUTs, we introduced compensation (Kc1-kc4 and Kc1.1-Kc4.1.; Equations (21)-(24); Section 3.2) to re-stabilise equilibrium and fitness and then generated new estimates and validated the model again ( Figure 4). . We considered the model to be accurate when the outcomes were included in the 95% CI. The second step was verification of this model by using data from Experiment 2 (exp 2, VIT E+) and Experiment 3 (exp 3, PUFA+). The VIT E+ scenario resulted in good accuracy (green box); thus, it did not require further adjustment. However, some OUTPUTs of the PUFA+ scenario were out of the 95% CI (red-coloured box) and showed unstable results. Hence, the model was improved through the introduction of compensatory values (Kc) obtained from intermediate OUTPUTs (Vit E/Vit Eox ratio). After compensation, the system was validated again (boxes surrounded by red dashed lines indicate > CI, whereas boxes surrounded by green dashed lines indicate < CI).
We set the values (in moles) of PUFA <55.5 and the ratio Vit E/Vit Eox <3.45 based on the results obtained during validation of the model (Section 3.1). We subsequently used and compared data from animal experiments (Section 2.3).

Experimental Plan: Animals and Dietary Treatments
Thirty blood samples were used to test and validate the model in two experimental sessions. Experiment 1 tested and verified the standard model (n = 10 samples), whereas validation was performed by using the same number of samples in Experiments 2 and 3 (n = 10 samples for each experiment) and diets with supra-nutritional levels of Vit E and PUFA [51]. Experiment 1.
Samples were collected from 10 New Zealand White male rabbits, 140 days old, reared at an experimental facility of Siena University and fed ad libitum a standard diet (Table A1 of Appendix A) containing 50 mg/kg alpha-tocopheryl acetate (vitamin-mineral premix) [52]. Experiments 2 and 3. We verified the model by using data from two distinct dietary interventions with high dosages of Vit E and PUFA: High VITAMIN E (VIT E+): 10 rabbits of the same age and breed mentioned above were fed ad libitum the standard diet (50 mg/kg alpha-tocopheryl acetate by vitaminmineral premix), including 200 mg/kg alpha-tocopheryl acetate. High PUFA (PUFA+): 10 rabbits of the same age and breed mentioned above were fed ad libitum a standard diet supplemented with 10% extruded flaxseed and 200 mg/kg alpha-tocopheryl acetate, following the rabbit requirements for a flaxseed enriched diet [53]. This diet had a 1.34-fold higher concentration of ALA than the standard diet.
The dietary protocols involved 50 days of adaptation during which the rabbits were only monitored and a subsequent 60 days during which blood samples were collected every 2 weeks.

Samples Collection
About 2 mL of blood was taken from the auricular marginal vein, after the local application of an anaesthetic cream (EMLA ® ), using a 2.5 mL syringe fitted with a butterfly needle. Serum was obtained from blood samples coagulated at room temperature for 2 h, whereas plasma was obtained from blood samples collected in tubes containing Na 2 -EDTA (ethylenediaminetetraacetic acid) and centrifuged immediately at 2500× g for 15 min at 4 • C.
The experimental trial was conducted in accordance with the Guiding Principles in the Use of Animals and approved by the Animal Ethics Committee of the Siena University (CEL AOUS; authorisation n • 265/2018-PR, ISOPRO 7DF19.23).

Analytical Determinations 2.5.1. Determination of the Fatty Acids
Lipids were extracted from the serum and feed according to the method of Folch et al. [54], and the esterification was carried out following the procedure of Christie [55].
The FA profile was determined by using a Varian gas chromatograph (CP-3800) equipped with a flame ionisation detector (FID) and a capillary column 100 m in length × 0.25 mm × 0.2 µm film (Supelco, Bellefonte, PA, USA). Helium was used as the carrier gas (0.6 mL/min). The split ratio was 1:20 for serum and 1:50 for feed. The oven temperature was programmed as reported by Mattioli et al. [56]. Individual fatty acid methyl esters (FAME) were identified by comparing the relative retention times of peaks in the sample with those of a standard mixture (FAME Mix Supelco; 4:0 to 24:0) plus cis-9, cis-12 C18:2; cis-9 cis-12 cis-15 C18:3; cis-9 cis-12 cis-15 C18:3 (all from Sigma-Aldrich, Steinheim am Albuch, Germany). The main PUFAs of both n-6 and n-3 series were quantified by using heneicosanoic acid methyl esters (Sigma Chemical Co., Schnelldorf, Germany) as an internal standard and are expressed in mg/mL of blood. Subsequently, each FA was converted to moles based on its specific molecular weight.

Determination of Vit E and AA
Vit E, in both the reduced form (α-tocopherol [α-TOH]) and oxidised form (αtocopherylquinone [α-TQ]), was quantified in plasma. Before lipid extraction, the samples were thawed at room temperature and homogenised by gentle pipetting.
α-TQ and α-TOH were quantified as described by Torquato et al. [20], with minor modifications. Briefly, 70 µL of plasma was spiked with α-TOH-d6 as an internal standard and mixed with 0.2 mL of sodium acetate buffer (pH 5) for the incubation (30 min at 37 • C) with 700 units of β-glucuronidase from Helix pomatia (Sigma-Aldrich, Steinheim am Albuch, Germany, G7017, Type HP-2). One hundred microlitres of glacial acetic acid and 1 mL of ethanol were added and, after mixing, the samples were incubated for 30 min at −20 • C and centrifugated at 3500 rpm for 5 min. Then, the ethanol was removed under a gentle stream of nitrogen and 2 mL of Milli-Q water was added to the sample. The lipid fraction was extracted twice by liquid/liquid partition with 3 mL of n-hexane/n-tert methyl tertiary-butyl ether (MTBE) mixture (2:1; v/v). After centrifugation at 3500 rpm for 5 min at 10 • C, the extracts were collected and dried under a stream of nitrogen. Then, the residues were resuspended with 80 µL of MOX reagent (Sigma Aldrich, TS-45950) and incubated at 70 • C for 2 h. The samples were further dried under nitrogen and then the analytes were derivatised with 50 µL of N-methyl-N-(trimethylsilyl)-trifluoroacetamide (Sigma-Aldrich, Steinheim am Albuch, Germany, 69479-5ML) in 50 µL of pyridine. Following incubation at 70 • C for 1 h, the samples were dried under a gentle stream of nitrogen and re-suspended in 125 µL of toluene for gas chromatography-mass spectrometry (GC-MS). Specifically, α-TOH and α-TQ were qualitatively and quantitatively assessed by using a mass spectrometer (Agilent Technologies, Milan, Italy; MSD 5975C, VLD with Triple Axis Detector) coupled to a gas chromatograph (GC 7890A) equipped with an Agilent VF-5 ms capillary column (15 m × 0.15 mm internal diameter [i.d.], 0.15 µm film thickness). The constant flow of the helium carrier gas was 1.0 mL/min and the injector operated in pulsed splitless mode (pressure 50 psi) at 300 • C. The oven temperature ramp for α-TOH and α-TQ was: from 120 • C (held for 1.5 min) to 240 • C at 113 • C/min, then to 270 • C at 8.5 • C/min (held for 1.5 min), to 350 • C at 28.3 • C/min, and finally the system was maintained under isothermal conditions for 5 min. The transfer line, ion source and quadrupole temperatures were set at 300, 230 and 150 • C, respectively. For each analyte and related recovery internal standard, the mass spectrometer operated in selected ion monitoring mode (EI+ SIM) at 70 eV: 508, 243 m/z for α-TOH-d6; 502, 237 and 430, 293 m/z, for α-TOH and α-TQ, respectively. We ensured high quality assurance and quality control of the performed methods by evaluating the testing linearity, detection limit, quantification limit, reproducibility and percentage of recovery values. The data are expressed as number of moles.
AA was measured in plasma by using a high-performance liquid chromatography system (HPLC) as described by Ross [57] with minor modifications as reported by Mattioli et al. [12]. AA was quantified by ultraviolet (UV) reverse-phase HPLC using a Waters 600 E System Controller HPLC equipped with a Waters Dual k 2487 detector (Milford, MA, USA) set at 262 nm. A 5 µm Ultrasphere ODS column (Beckman, San Ramon, CA, USA) was used with acetonitrile:water (49:51, v/v) as the mobile phase at a flow rate of 0.8 mL/min. AA concentrations were calculated by peak areas determined using an Agilent 3395 integrator (Agilent Technologies, Santa Clara, CA, USA); the results are expressed in moles.

Determination of ROS, Malondialdehyde, F 2 -Isoprostane, F 3 -Isoprostane, and F 4 -Neuroprostane
ROS levels in plasma were evaluated by using a commercial kit (Diacron, Grosseto, Italy) and are expressed as moles of H 2 O 2 .
Lipid peroxidation in plasma was assessed by the MDA level using malondialdehyde tetrabutylammonium as a standard (Sigma-Aldrich, Steinheim am Albuch, Germany) and indicated in the model as Aldehydes. Rabbit plasma samples were mixed in a 0.04 M K + -phosphate buffer (pH 7.4) containing 0.01% BHT (1:5 w/v, 0 • C) to prevent artificial oxidation of free PUFA during the assay. The supernatants resulting from deproteinisation with acetonitrile (1:1) were used for aldehyde analysis after pre-column derivatisation with 2,4-dinitrophenylhydrazine [58]. The samples were stirred immediately, extracted with 5 mL of pentane and dried under a nitrogen flow. The aldehyde hydrazone was quantified by isocratic HPLC using a Waters 600 E system controller HPLC instrument equipped with a Waters Dual λ 2487 UV detector set at 307 nm. A 5 µm Ultrasphere ODS C18 column was used with a mobile phase composed of acetonitrile (45%) and 0.01 N HCl (55%) at a flow rate of 0.8 mL/min. The aldehyde concentration was calculated based on the peak areas using an Agilent 3395 integrator. The results are expressed as moles of Aldehyde.
The concentration of all the blood traits was divided by 10 16 to make the data easier to read.

Statistical Evaluation
A linear model (SAS, [60]) was used to evaluate the fixed effect of diet (Standard, VIT E+, PUFA+) on blood traits resulting from chemical evaluation. We used multiple comparisons (Bonferrini's test) to establish the significance of differences between dietary treatments (p < 0.05). We calculated the mean, standard error and 95% CI of the analytical data.
VENISM software (release 9.1.0; Ventana Systems, UK Ltd., Wiltshire, UK) created several non-linear equations for estimating the outcomes at different lag-times based on the relationship between variables, which were imposed by the developer. We evaluated the accuracy of the models by comparing the 95% CI of the analytical data with the relative outputs estimated for calibration and then for the validation procedure.
As exemplified in the cognitive map (Figure 4), for Experiment 3 when the outcomes were far from the CI (e.g., PUFA+; Table 1), the model became unstable and needed compensation. We calculated this compensation from the effect of PUFA increase on the Vit E/Vit Eox ratio by fitting a polynomial regression curve. We used a first derivative equal to 0 to calculate the flex points (average, minimum and maximum) and we used the relative y value to assess the compensation rates (Kc1-Kc4 and Kc1.1-Kc4.1) (Equations (21)-(24)). We re-evaluated the accuracy of the resulting model by using the above-mentioned procedure.

Model Validation in the VIT E+ and PUFA+ Scenarios
We compared the analytical and model outcomes for the VIT E+ and PUFA+ scenarios. As shown in the cognitive map (Figure 4), the VIT E+ scenario was estimated accurately with the standard model, whereas the PUFA+ scenario results were far from the 95% CI.
The VIT E+ scenario resulted in similar amounts of IsoNeuroPs and AA and lower Aldehydes than the standard one. The Vit E concentration resulting from addition to the diets was similar to what was recorded in the standard scenario (458.41 × 10 −5 vs. 422.00 × 10 −5 moles, respectively), but its oxidised form (Vit Eox) was 3-fold higher.
On the contrary, when PUFA INPUTs increased, the model outcomes became unstable ( Figure 6) and far from the 95% CI (data not shown; IsoNeuroPs = 0.34 moles, Aldehydes = 16 moles, Vit Eox = 32.2 × 10 −5 moles). In this situation, one of the most relevant disturbing factors was the Vit E/Vit Eox ratio, which decreased as PUFA increased. The flex points of this polynomial relationship correspond to about 55.5 moles PUFA and a Vit E/Vit Eox ratio of 3.45 ( Figure 7). As mentioned previously, we applied some compensatory rates (for Kc1.1-Kc1.4, 12.50, 0.10, 37.4 and 0.50, respectively) in this scenario to render the molecular flows more physiological and the results more accurate. After compensation, we re-validated the system (Table 1), and almost all the estimated values (except for IsoNeuroPs and Vit E) were within the 95% CI. This scenario showed about 13-fold more PUFA than the others, and it also showed higher ROS, whereas Vit E was quite similar. Supra-nutritional PUFA resulted in an almost 3-fold higher level of IsoNeuroPs and Aldehydes compared with the standard model.

Discussion
Organisms exhibit a myriad of homeostatic mechanisms based on internal and external interactions, which modulate biological responses at different scales and dynamical rates [61]. Building a mathematical model can be an effective way to describe and understand some aspects of such variability. The best outcome of such models is to predict the behaviour of these complex systems, returning values close to the measured one. The rationale of modelling is based on the simplification of the system by assuming some approximations and uncertainty. Moreover, to better adhere to in vivo values, the model should approximate the complexity and the number of variables of the biological system as much as possible. According to such methodological needs, in the present study we proposed a modelling strategy for a living system, using as a case study the metabolism of fat-soluble nutrients in the rabbit. Within this framework, the dynamic model outlined some crucial aspects in the balance of the ROS-PUFA-antioxidant system. In this model system, we evaluated some common OUTPUTs that approximate the oxidative status of the investigated organism, although the quantification of some of them (i.e., ROS, aldehydes, isoprostanoids) is not fully representative of oxidative stress and requires mentioning their specific chemical entities [62]. Starting from a standard condition (standard diet), we introduced dietary changes in the model system (VIT E+ and PUFA+ scenarios) and predicted the relative OUTPUTs (ROS, PUFA, Vit E, AA, IsoNeu-roPs, Aldehydes and Vit Eox). We estimated the accuracy of the prediction by comparing analytical data obtained in the corresponding experimental models of nutritional intervention ( Figure 4, Table 1). The accuracy of the model was high and the OUTPUT variables were well described by the simulation in both the standard and VIT E+ scenarios, but not in the PUFA+ scenario.
VIT E+ resulted in IsoNeuroP and AA levels similar to the standard model; on the contrary, the Aldehyde levels were reduced. Vit E levels were increased slightly, but its oxidised form (Vit Eox) was 3-fold higher in this comparison. As expected, the addition of 200 mg/kg of α-tocopheryl acetate (VIT E+) reduced the global oxidation index (Aldehydes) of the standard model without affecting specific oxidation markers (IsoNeuroPs). Although excess Vit E may lead to undesirable effects, such as accumulation of this vitamin in fat depots [20] and its increased exposure to autoxidation processes [24,29], based on the formation of α-TQ, the system can utilise a certain amount of the administered Vit E to donate H atoms to other substrates. This finding confirms its redox interaction with the other lipids of the system. In fact, the Vit E levels remained in the physiological range [24] in VIT E+ conditions, that is, similar to the levels observed in animals fed a standard diet (458.41 × 10 −5 vs. 422.00 × 10 −5 moles, respectively).
The PUFA+ scenario presented a marked increase in blood PUFA (about 13 times higher than the other scenarios). This elevation apparently destabilised the main relationship between key molecules involved in detoxification of PUFA-derived oxidative products (i.e., Vit E and its oxidised form). Accordingly, when PUFA INPUT was >55.5 moles (e.g., 80 moles), the model became unstable (Figure 6) and the outcomes lacked accuracy. Hence, we introduced compensatory rates in the model. This type of correction is not a mathematical artifact due to fitting analytical and estimated values, but rather a way to keep considering metabolic changes that occur in a system that becomes unstable.
It is widely known that high dietary levels of PUFA increase the risk of autoxidation and could impair the antioxidant defence of the organism [12,63]. Moreover, deficient antioxidant defence in a tissue could result in higher risk of chronic inflammation and development of degenerative diseases [64]. Obviously, increased PUFA intake was well tolerated in the healthy rabbits utilised as the experimental model in this study. This elevated intake stimulated compensatory responses relevant to the redox homeostasis of the animal. Considering these facts, we introduced different compensatory rates (Kc1-Kc4 in the PUFA+ scenario. These compensatory rates represent critical checkpoints of the system; when self-adjustments are required (e.g., change of flow or genetic reprogramming), this indicates unbalanced homeostasis and the need for appropriate corrections to be introduced in the system (i.e., changes in Vit E/PUFA intake or the reduced glutathione [GSH]/oxidised glutathione [GSSG] ratio, etc.). The reasons for each compensatory factor manually introduced in the system are as follows: Kc1 aims to counteract isoprostanoid accumulation that is activated when PUFAs concentrations exceed physiological levels. Indeed, in healthy conditions, isoprostanoids circulate in the plasma and are excreted in the urine [27]. However, when these metabolites reach certain levels, they undergo rearrangement, dehydration or conjugation with glucuronide in the liver to yield a variety of secondary metabolites, thus leading to a reduction in their concentrations and the risk of adverse effects [6].
Kc3 intervenes when the levels of the free radical-derived metabolites of Vit E, α-TQ or the same Vit E, are too high (Figure 6), implying increased expression of detoxifying genes encoding CYP450s involved in the metabolism of this vitamin and other antioxidants, as well as of xenobiotics, electrophiles, heavy metals, etc. [66].
Finally, Kc4 compensates for external perturbation variables (i.e., diet), eventually interfering with antioxidant capacity and with oxidative processes.
Independently of these compensations, the PUFA+ scenario showed more PUFA and ROS, whereas Vit E was quite similar to the other models (p < 0.05). However, IsoNeu-roPs, Aldehydes and Vit E ox were almost 2-3-fold higher compared with the standard model (Table 1). Although the introduction of these compensatory rates improved the accuracy of the PUFA+ scenario, some traits (Vit E, IsoNeuroPs) remained slightly outside the 95% CI. The differences observed could be ascribed to different factors: on the one hand, the sensitivity of analytical methods, including systematic errors of the measurements, are not considered in the model; on the other hand, the model has limits because it forcefully misses or miscalculates some molecular interactions and their outcomes. Furthermore, the accuracy of the model should improve when the number of samples increases.
In this view, comparison of in vivo data with what is estimated by the model highlights some advantages and disadvantages of the modelling strategy. Regarding advantages, the model provides an integrated tool to assess the time-course of a complex set of variables representing a comprehensive and dynamic picture of the system. This information cannot be obtained from punctual analyses and techniques that describe the system only with static snapshots. Regarding disadvantages, the model is naturally based on limitations and assumptions dealing with the complexity of the aspects considered in the study. Furthermore, some of them are still unknown or undetermined.
There are only a few examples in the literature of similar procedures for modelling biological processes. Musesti et al. [67] and Giantesio et al. [68] created a mathematical model of muscle tissue modification to study cell ageing in sarcopenia, using physiological parameters (i.e., muscle mass, electrical impedance) and variations in external factors (e.g., diet, physical activity, pharmacological treatments, environmental pollution exposure, etc.). They applied the model to analyse two typical symptoms of sarcopenia [68]: loss of mass and loss of force, finding non-linear functions that described these alterations of muscle physiology. These authors also applied a compensatory strategy. Specifically, they developed these functions by defining some constants (identified in their studies as d) representing the constant damage associated with the different symptoms of sarcopenia. This approach demonstrated the relevance of modulating these constants in the simulation to depict different scenarios of the biological system.
Although our approach may have biological and mathematical limits, it is apparent that it could be very useful when applied to determine single or combined effects of dietary interventions on the redox homeostasis of the rabbit. In particular, the scientific and technical literature has shown that rabbit bucks and/or does require high levels of PUFA and Vit E for suitable reproductive performance [53,69]. Indeed, germ cells have a particularly high level of LC-PUFA in their membranes. Thus, to ensure good fertility, it is necessary to include certain amounts of PUFA (i.e., LA, ALA, EPA, DPA, and DHA) and antioxidants (i.e., Vit E and AA) in their diets [14]. However, the levels of these compounds should be regulated objectively and carefully; otherwise, they could miss or even hinder the expected improvements. For these reasons, the development of modelling tools holds great potential for precision nutrition studies using animal models and humans.

Conclusions
We developed a mathematical model based on computer software that mimics the complex interactions of the ROS-PUFA-antioxidant system in rabbit, namely the VENSIM model. This can be used to estimate a priori the effect of dietary treatments (i.e., antioxidants, PUFAs) on the redox homeostasis of the rabbit, and in turn on the animal's health.
In this view, the model can assist in assessing the effect of oxidative challenges and dietary interventions with bioactive lipids such as Vit E and PUFA on the reproductive efficiency of animals/humans, which is a physiological function very sensitive to redox imbalance and oxidative stress. This model system can be implemented and easily be instructed to simulate a wide panel of processes and to assess several phenotypic traits.
In particular, self-learning processes able to identify potential health risks linked to imbalances of the system can be developed. Further improvements of the model would permit studying other systems (i.e., male/female reproduction, key organs, etc.) and the response to exogenous factors different from diets, such as environmental, pathological and pharmacological perturbations.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Formulation, proximate analysis and fatty acid profile of the control (Standard) and polyunsaturated + (PUFA+) and vitamin E + diet (Vit E+). Fatty acids provided by flaxseed C18:3 ALA g/kg 9.23 9.05 12.08