Dynamic simulation of cardiolipin remodeling: greasing the wheels for an interpretative approach to lipidomics.

Cardiolipin is a class of mitochondrial specific phospholipid, which is intricately involved in mitochondrial functionality. Differences in cardiolipin species exist in a variety of tissues and diseases. It has been demonstrated that the cardiolipin profile is a key modulator of the functions of many mitochondrial proteins. However, the chemical mechanism(s) leading to normal and/or pathological distribution of cardiolipin species remain elusive. Herein, we describe a novel approach for investigating the molecular mechanism of cardiolipin remodeling through a dynamic simulation. This approach applied data from shotgun lipidomic analyses of the heart, liver, brain, and lung mitochondrial lipidomes to model cardiolipin remodeling, including relative content, regiospecificity, and isomeric composition of cardiolipin species. Generated cardiolipin profiles were nearly identical to those determined by shotgun lipidomics. Importantly, the simulated isomeric compositions of cardiolipin species were further substantiated through product ion analysis. Finally, unique enzymatic activities involved in cardiolipin remodeling were assessed from the parameters used in the dynamic simulation of cardiolipin profiles. Collectively, we described, verified, and demonstrated a novel approach by integrating both lipidomic analysis and dynamic simulation to study cardiolipin biology. We believe this study provides a foundation to investigate cardiolipin metabolism and bioenergetic homeostasis in normal and disease states.

concentration of isomeric cardiolipin species, which have been identifi ed by MDMS-SL, but also predicts the existence of a variety of very low abundance cardiolipin molecular species, many of which cannot be accurately determined using currently available technologies. Furthermore, experimentally measured profi les of cardiolipin molecular species can be recapitulated with this simulation utilizing specifi c distributions of acyl chains refl ecting the coordinated activities of phospholipases, acyltransferase(s), and/or transacylase(s). Collectively, we demonstrated that the cardiolipin dynamic simulation combined with MDMS-SL analyses of PG, phosphatidylcholine (PC), phosphatidylethanolamine (PE), and acyl CoA is a powerful new approach for determining the presence of isomeric cardiolipin molecular species and assessing the enzymatic activities that meditate the extent of cardiolipin remodeling. We contend that this dynamic simulation will provide novel insight into the mechanisms regulating the complex tissue-specifi c cardiolipin molecular species distributions and that new insight into the molecular mechanisms underlying alterations in acyl chain selectivity in health and disease can be accrued.

Mitochondrial isolation
C57BL/6J wild-type male mice (4 months of age) were purchased from the National Institute on Aging (Bethesda, MD). All animal procedures were performed in accordance with the Guide for the Care and Use of Laboratory Animals and were approved by the Animal Studies Committee at Washington University. Mice were euthanized by asphyxiation with carbon dioxide. Mouse organs were perfused with phosphate-buffered saline plus heparin to remove possible plasma contamination before dissection. Brain cortex, heart, liver, and lung of mice were harvested immediately after perfusion. Nonsynaptic brain mitochondria were isolated as previously described in detail using discontinuous Ficoll and sucrose gradients ( 25,26 ). Myocardial, hepatic, and lung mitochondria were isolated as previously described ( 27 ). All experiments were performed on ice in a 4°C cold room. Briefl y, the heart, liver, or lung was diced using a razor blade and homogenized in homogenization buffer [HM; 0.21 M mannitol, 70 mM sucrose, 0.1 mM EDTA-K, 0.5% BSA, 10 mM Tris-HCl, 1× protease inhibitor cocktail (complete®, Roche, Germany) pH 7.4] in a ratiometric volume of HM buffer to tissue of 8:1. The homogenate was then centrifuged at 700 g for 5 min. The supernatant was collected and the pellet was homogenized and recentrifuged. The supernatant was collected, and the pooled supernatants were centrifuged at 900 g for 5 min. The pellet was discarded, and the supernatant was collected and centrifuged at 10,000 g for 10 min. The pellet was resuspended in HM buffer and layered on a 1.6, 1.2, 1.0, and 0.8 M discontinuous sucrose gradient and spun for 50,000 g for 2 h in a SW 28 rotor. The highly enriched mitochondrial fraction was collected at the 1.6/1.2 M sucrose interface and recentrifuged at 18,000 g for 15 min. The pellet was resuspended in HM buffer without BSA and recentrifuged at 8,000 g for 10 min. The pellet was resuspended in HM buffer without BSA, and the protein concentration of the resuspended mitochondrial fractions was determined using the BCA protein assay (Thermo Scientifi c, Rockford, IL). disease states ( 12,13,(15)(16)(17)(18)(19). This diversity is regulated by differences in cardiolipin biosynthesis and catabolism, remodeling acyl selectivities, and transacylase as well as acyltransferase activities ( 20,21 ). Cardiolipin is de novo synthesized in the mitochondrial membrane by the condensation of phosphatidylglycerol (PG) and cytidine diphosphate-diacylglcyerol (CDP-DAG) ( Fig. 1 ). The newly formed polyglycerophospholipid is referred to as immature cardiolipin, which is then deacylated and extensively remodeled using acyltransferase(s) and/or transacylase(s) ( 21,22 ). Cardiolipin remodeling utilizes acyl CoA and acyl chains from the sn -2 position of choline glycerophospholipids (PC ) and ethanolamine glycerophospholipids (PE) accordingly as acyl donors ( 20,23,24 ). The balance of acyltransferase and transacylase activities as well as acyl chain selectivities is refl ected in the acyl chain pool, which sculpts nascent cardiolipins into mature cardiolipin molecular species.
Based on known biological pathways and MDMS-SLderived data, we developed a dynamic simulation approach that recapitulates the temporal reconstruction of immature cardiolipins into mature cardiolipins by remodeling from available donor acyl chains. The sequential remodeling events can be visualized graphically for intuitive comparison of the dynamic aspects of cardiolipin remodeling. The biological parameters in these simulations can be modifi ed to compare cardiolipin profi les obtained from MDMS-SL. Moreover, the algorithmic approach utilized in this simulation not only demonstrates the presence and The schematic cardiolipin biosynthesis and remodeling pathway. CDP-DAG is used in the mitochondria to form the metabolic intermediate PG phosphate for the synthesis of PG through a dephosphorylation activity. Newly synthesized cardiolipin (immature cardiolipin) is formed by the condensation of PG and CDP-DAG. Immature cardiolipin is then deacylated to form monolysocardiolipin and then reacylated using acyl chains from acyl CoA or the sn -2 acyl chain of phosphatidylcholine and phosphatidylethanolamine. This remodeling process yields mature cardiolipin. producibility of 90% and a dynamic range of 100-fold were achieved by using this method.

Mathematical computational strategy for simulating cardiolipin remodeling
Generation of immature cardiolipin molecular species. Cardiolipin is de novo synthesized by the condensation of PG and CDP-DAG. CDP-DAG is utilized for the synthesis of PG in the mitochondria ( Fig. 1 ) and thus should contain similar molecular species to PG profi les, which has been previously demonstrated ( 37 ). Due to diffi culties in analyzing mitochondrial CDP-DAG, only PG molecular species were utilized to generate immature cardiolipin molecular species. PG molecular species were quantifi ed, and the sn -1 and sn -2 acyl chain composition was determined by MDMS-SL analysis ( 5 ). Based on the composition of PG species, newly synthesized cardiolipin molecular species were generated using a customized PERL script. The sn -1, sn -2, sn -1', and sn -2' positions of immature cardiolipin species were formed based on the sn -1 and sn -2 position and concentrations of PG molecular species ( Fig. 2 ) ( 1 ). The generated distribution provided both the initial condition for the dynamic simulation as well as the distribution of immature cardiolipin that would be incorporated in the system at a constant rate.
Modeling the rate, catabolism, and sn-1 and sn-2 selective remodeling of cardiolipin molecular species. To quantitatively describe the cardiolipin remodeling process, we developed software that simulates the dynamics of transferring acyl chains into cardiolipin. The simulation is based on incorporation of determined fatty acids from sn-2 positions of PC and PE or acyl CoA at designated concentrations into cardiolipin molecular species under mass action kinetics, with all four cardiolipin acyl chain positions treated equivalently. We fi rst describe the rate of change of a given cardiolipin species with the chains ␣ 1 ␣ 2 ␣ 3 ␣ 4 at the sn -1, sn -2, sn -1', and sn -2' positions, respectively. The process of chain substitution from PC, PE, and acyl CoA into cardiolipin corresponds to a dynamic equation : Here ( 1) sn PC k is the reaction rate constant between PC and one of the sn -1 positions of cardiolipin, and ( 2) sn PC k is the reaction rate constant between PC and one of the sn -2 positions. Thus, the model can accommodate differential kinetic fl ux present at the sn -1 and sn -1' positions versus sn -2 and sn -2' positions in cardiolipin that mimic the phospholipase activities that contribute to the generation of monolysocardiolipin and cardiolipin remodeling. For the analogous PE terms, the PE-specifi c values ( 1) sn PE k and ( 2) sn PE k replace the corresponding PC rate parameters and likewise for the acyl CoA terms. For simplicity, our model assumes that the relative reactivity of the sn -1 and sn -2 positions is a property of the cardiolipin molecule, i.e., ( 2) ( 1)

Instrumentation and MS
A triple-quadrupole mass spectrometer (Thermo Scientifi c TSQ Quantum Ultra Plus, San Jose, CA) equipped with a Nanomate device (Advion Biosciences Ltd, Ithaca, NY) and Xcalibur system software was utilized for MDMS-SL analysis as previously described ( 28,29 ). Lipid molecular species were identifi ed, quantifi ed, and acyl chain designation determined using a MDMS-SL approach as previously described ( 1,5,30,31 ). Product-ion analyses of cardiolipin molecular species by using a high mass accuracy/high resolution mass spectrometer [i.e., Orbitrap-LTQ mass spectrometer (Thermo Scientifi c San Jose, CA)] were also preformed as previously described ( 32 ) to verify the presence of very low abundance cardiolipin isomers.

Acyl CoA extraction and analysis
Extraction and analysis of acyl CoA was performed as previously described ( 34 ) with modifi cation. Briefl y, 4-month-old C57BL/6J male mice were sacrifi ced, and tissue samples were immediately harvested and clamp frozen in liquid nitrogen. The rapid rate of tissue harvest ( ‫ف‬ 15-30 s) is critical for accurate assessment of acyl CoA levels and distribution. Tissue samples were then pulverized in liquid nitrogen, and 50-80 mg of tissue powder was weighed and homogenized in chloroform/methanol (1:2, v/v) in an ice/ethyl alcohol bath to immediately denature enzymes. An internal standard was added (17:0 CoA) at a concentration of 50, 90, 50, and 60 pmol/mg protein for the liver, heart, lung, and cortex, respectively. This internal standard was selected, because the endogenous mass of this species was minimal, as demonstrated by MALDI-time of fl ight (TOF)/MS lipid analysis (supplementary Fig. I). Chloroform and deionized water were added, yielding a ratio of chloroform/methanol/deionized water of 2:2:1 (v/v/v) and a volume of at least 10 l/mg of tissue powder. The samples were then homogenized in an ice/ethyl alcohol bath for 2 min and then incubated on the ice/ethyl alcohol bath for 30 min, vortexing occasionally. The homogenate was then centrifuged at 23,000 g for 10 min in a 4°C cold room. The methanol and deionized water layer (upper phase) was collected. Remaining tissue residue was reextracted using the same procedure. The combined aqueous phases were dried under nitrogen and reconstituted in 50 l of deionized water. A mixture of 5 l of acyl CoA extract and 5 l of 9-aminoacridine solution [10 mg of 9-AA/ml in acetone/methanol (98:2, v/v)] was spotted on a MALDI plate and analyzed using a MALDI-TOF mass spectrometer ( 34 ). Quantitative analyses of acyl CoA species by MALDI-TOF-MS with 9-aminoacridine as matrix were similarly assessed as previously described ( 35,36 ) after 13 C-deisotoping ( 1 ). A re-tions we use a timestep of ⌬ t ‫ف‬ 0.05-0.15. Whenever practical, a larger value of ⌬ t is used to reduce the simulation computation time.
The dynamic simulation can determine the cardiolipin concentrations at each time point in three output formats: 1) the concentrations of all molecular species as specifi ed by the acyl chains in each of the four acyl positions (full chain output), 2) the aggregate concentrations of all molecular species at each designated number of carbons and double bonds, a format that can be easily compared with experimental mass spectrometry measurements (mass spectrometry output), or 3) a hybrid format in which the cardiolipin species are sorted by the numbers of carbons and double bonds, which allows for simplifi ed analysis of the isomeric species that correspond to each ion of cardiolipin species determined at a certain m/z by MS.
Acyl chain distributions used for cardiolipin remodeling. The acyl chains utilized for the formation of mature cardiolipin molecular species could originate from a transacylation reaction from PC or PE or via an acyltransferase reaction from acyl CoA (38)(39)(40). To model these selective aspects of the remodeling machinery, we devised an approach that allows for the input of PC and PE sn -2 positional acyl chain composition as well as acyl CoA chain data utilizing designated acyl chain selectivity ( Fig. 2) ( 3 ). Alternatively, it is possible to manually select the acyl chains from PC and PE as well as from acyl CoA to be contributed to cardiolipin remodeling. This aspect of the remodeling program allows for the capability of defi ning the possible ratios of acyltransferase versus transacylase activities based on the fatty acids that are present in discrete cardiolipin molecular species, in addition to the interpretation of acyl selectivity.
To account for the steady synthesis of immature cardiolipin from PG, which will replace the proportion of cardiolipin that is catabolized, we added two additional terms to the simulation. One term corresponds to the constant catabolism of a small fraction of the cardiolipin, described by an exponential decay parameter ␥ . This catabolism is assumed to be nonspecifi c, i.e., all cardiolipin species are catabolized at the same rate. The second term corresponds to a constant injection of immature cardiolipin that keeps the total quantity of cardiolipin in the system constant. Inclusion of these terms allows for the simulation of the relative importance of rates of cardiolipin turnover (catabolism and de novo synthesis) versus rates of acyl chain substitution. Overall, the equation describing the kinetics is given by: where the chain exchange terms are described in the equation above. This equation is integrated numerically with time step ⌬ t , a parameter that determines the amount of cardiolipin that will be remodeled in a single iteration using designated acyl chains from PC and PE sn -2 positions, as well as acyl CoA. In practice, the numerical integration is performed using dimensionless variables, with a timestep small enough that only a small fraction of the cardiolipin is remodeled in any given iteration to maintain realistic behavior. In these dimensionless variables, this condition is typically satisfi ed as long as ⌬ t р 1 (with some dependence on the distribution of species concentrations of each of the acyl chain sources), and as a conservative criterion in most simula-

Fig. 2.
The strategy employed for cardiolipin dynamic simulation. 1: PG molecular species including regiospecifi city were identifi ed and quantifi ed. Predicted immature cardiolipin was generated based on the concentration of PG species. 2: Immature cardiolipin was then deacylated with all chains being equally selected or with sn -1 or sn -2 selectivity followed by the reacylation with acyl chains 3: originating from phosphatidylcholine and phosphatidylethanolamine sn -2 positions (a) or acyl CoA (b ) (with acyl selectivity inferred) representing transacylase and acyltransferase activities, respectively. Additionally, at this point, a portion of the cardiolipin molecular species can be nonselectively catabolized and immature cardiolipin species can be constantly added to represent a proportion of cardiolipin biosynthesis. 4: Cardiolipin remodeling will reach an equilibrium state in which the distribution of cardiolipin molecular species will not dramatically change relative to the distribution of acyl chains used for remodeling. cardial mitochondrial PG (supplementary Table II) revealed that there is only a minimal contribution of cardiolipin molecular species originating from de novo biosynthesis in the mature remodeled myocardial cardiolipin distributions identifi ed by MDMS-SL ( Fig. 3A ). Thus, the amount of decay and replacement estimated in the dynamic simulation for myocardial cardiolipin was minimal due to the lack of immature cardiolipin species in the mature cardiolipin distribution. We determined by MDMS-SL as well as experimentation a distribution of acyl chains as 16:0 (0.27%), 16:1 (0.6%), 16:2 (0.03%), 18:0 (0.13%), 18:1 (9.92%), 18:2 (67.79%), 18:3 (0.19%), 20:1 (0.09%), 20:2 (1.98%), 20:3 (2.1%), 20:4 (0.84%), 22:5 (1.43%), and 22:6 (14.36%) and used this distribution to simulate myocardial cardiolipin remodeling that generated a cardiolipin profi le best representing the profi le obtained by MDMS-SL. A ⌬ t value (see "Materials and Methods") of 0.1 was used in the simulation. Figure 3 demonstrates the sequential remodeling of cardiolipin from the predicted immature cardiolipin distribution after 0, 50, 100, and 250 iterations ( Fig. 3A-D , respectively). Simulation data can be found in supplementary materials and viewed by downloading the customized cardiolipin remodeling viewer.
During the progression of remodeling steps, the stepwise iterations demonstrate the structural remodeling of cardiolipin molecular species at different stages of the remodeling process. The process of cardiolipin remodeling from predicted immature state to a mature state involves catabolism of the majority of species in immature cardiolipin ( Fig. 3A ), which then progress into a broad distribution of cardiolipin molecular species ( Fig. 3B ). For example, a similar broad distribution of cardiolipin molecular species was observed between day 5 and 15 for normal mouse cardiac development (data not shown) as well as in mature mouse brain ( 16,26 ). It is interesting that this broad distribution of cardiolipin molecular species is also observed in various disease states, culture conditions, and transgenic models of enzymes that catalyze cardiolipin remodeling ( 23,(41)(42)(43)(44)(45). This suggests a common remodeling event or role of various acyltransferases, transacylases, or phospholipases involved in this stage of the remodeling process. After 100 iterations of in silico remodeling, the mature myocardial cardiolipin distribution becomes more evident ( Fig. 3C ) and reaches an equilibrium state after 250 iterations ( Fig. 3D ). This molecular species distribution corresponds closely to the distribution of cardiolipin molecular species obtained from MDMS-SL in mouse myocardial mitochondria ( Fig. 3E , Table 1 ). Similarities in designated molecular species can also be viewed in the hybrid output of the dynamic simulation for the heart (supplementary materials), which demonstrates the accuracy in identification of the molecular species obtained by MDMS-SL and the dynamic simulation approach. In all, the dynamic simulation output demonstrates that the level of detection of cardiolipin molecular species by current technologies is around the threshold of 0.1% of total cardiolipin from myocardial mitochondrial extracts. The dynamic simulation suggests other extremely low abundance species that likely exist in the cardiolipin distribution but are below the detection limit of MDMS-SL under the experimental conditions ( Table 1 ).
Visualization and dynamic simulation output for cardiolipin remodeling. To view the sequential remodeling of cardiolipin molecular species, a graphical viewer was written in JAVA for enhanced portability across operating systems. This viewer displays the percentage distribution of cardiolipin species as organized by the carbon:double bond composition (an output of mass spectrometric format). Sequential remodeling of cardiolipin species can be viewed on a frame-by-frame basis or as an animated movie. All software is available at http://bioinformatics.bc.edu/~jchuang/dynamic.html.

RESULTS
Analysis of cardiolipin molecular species and their acyl chain regioisomers in myocardial, hepatic, lung, and brain mitochondria revealed a diverse distribution that provided insight into the biochemical processes that mediated the distinct distributions of cardiolipins in these tissues (supplementary Tables I and VII). Based on the premise that the diversity of cardiolipin molecular species in various organs was dependent on the availability of different acyl chains in PC and PE as well as the acyl CoA composition and the acyl selectivities of the involved acyltransferases and transacylases, a dynamic model was constructed to identify the critical reactions necessary for the generation of mature cardiolipins.
Using the MDMS-SL measured concentrations of individual cardiolipin species as well as those of the donor acyl chains from PC, PE, and acyl CoA, we computed the relative contributions from each of the donor substrates to the observed cardiolipin distribution. A distribution of acyl chains was generated from PC, PE, and acyl CoA that resembled the acyl composition of myocardial, hepatic, lung, or brain cardiolipin, and this distribution was used for the chain substitution process (see "Materials and Methods"). Based on the analysis of predicted immature cardiolipin from PG molecular species, the proportional amount of newly synthesized cardiolipin was inferred to contribute to the observed cardiolipin distribution in each tissue. Generally, the amount as well as ratio of saturated acyl chains found in PG molecular species served as reliable indicators of the relative content of immature cardiolipin molecular species being represented in the mature pattern of remodeled cardiolipin. This assumption seems reasonably justifi ed, because the content of saturated acyl chains in mature cardiolipins are quite low and these are likely derived from either de novo cardiolipin biosynthesis or from saturated acyl chains derived from acyl CoA during remodeling. Because the sn -2 positions of PC and PE contain almost solely unsaturated acyl chain species, it is unlikely that the saturated acyl chains in mature cardiolipin transacylated from either the PC or PE pools. Using this approach, we examined the dynamic time course of altered cardiolipin acyl chain distributions of four bioenergetically diverse organs: the heart, liver, and brain (see next sections), as well as lung (see supplementary materials).

Dynamic simulation of myocardial cardiolipin molecular species
Examination of the predicted immature cardiolipin molecular species generated by the condensation of myo-  (15). The mouse myocardial cardiolipin molecular species detected and quantifi ed by MDMS-SL as well as determined by the dynamic simulation were tabulated in Table 1 . Mouse myocardial mitochondrial lipids were extracted by a modifi ed Bligh and Dyer procedure as described in the section of "Materials and Methods." The cardiolipin molecular species in the lipid extracts were identifi ed and quantifi ed as previously described ( 15 ). The mass content and distribution of cardiolipin species represent the mean of four separately prepared mitochondria from different animals. The distribution of cardiolipin molecular species determined by the dynamic simulation represents the one frame output 0.1 ␣ 250 iterations for heart cardiolipin. This data was obtained using the MS output function. Values < 0.01 percent are not displayed but can be found in supplementary tables. CL, cardiolipin.

Dynamic simulation of brain cardiolipin molecular species
The analysis of the distribution of brain cardiolipin molecular species is more complicated than that of myocardial or hepatic cardiolipin. This is due to the broad distribution, the complex acyl chain composition, and the presence of numerous isomeric molecular species compared with the cardiolipin species present in the other organs. Additionally, immature cardiolipin molecular species formed from PG could account for a portion, small but signifi cant, relative to the mature cardiolipin distribution ( Fig. 5 ). Thus, the residual portion of cardiolipin biosynthesis being represented in the overall pattern would complicate the analysis of predicted species being generated in the overall distribution, depending on the rate and stage of remodeling.
To incorporate the proportion of immature cardiolipin species from de novo cardiolipin biosynthesis into the cardiolipin remodeling simulation, a decay rate of 0.01 was used with an initial acyl chain distribution of 16 ). This distribution generated the best-fi t cardiolipin distribution compared with that obtained by MDMS-SL. A ⌬ t of 0.1 was used in the simulation. The sequential remodeling of brain cardiolipin using these designated acyl chains led to a greater abundance of species that contain longer chains and greater unsaturation ( Fig. 5B, C ) compared with the distribution of heart and liver cardiolipin. Experimental attempts to simulate brain cardiolipin molecular species with a lower decay rate than 0.01 resulted in less diversity of molecular species and an increase in the percentage distribution of the major molecular species (data not shown) that did not accordingly fi t the actual mature cardiolipin distribution determined by MDMS-SL ( Fig. 5D, E ). Again, the profi le generated by the dynamic simulation resembled the percentage distribution determined by MDMS-SL ( Table 3 ) and maintained an approximate threshold of detection around 0.15% of total brain cardiolipin.

Validation of isomeric cardiolipin species determined by the dynamic simulation
To substantiate the accuracy of our dynamic simulation approach, we performed product ion analysis of ions containing very low abundance cardiolipin isomers by using a high mass accuracy/resolution mass spectrometer, which could not be detected by MDMS-SL at its current stage but were determined by our simulation. Owing to the low background noise, high sensitivity, and high mass accuracy/resolution, the presence of very low abundance isomers of an ion can be identifi ed through high mass accuracy/resolution mass spectrometry as demonstrated previously ( 32 ).
For example, selection and collision-induced dissociation of the cardiolipin species at m/z 723.48 present in mouse myocardial cardiolipins demonstrated an abundant fragment at m/z 279.23 in the fatty acyl fragment The equilibrium reached in the simulation is dynamic, in which remodeling maintains the cardiolipin distribution, and acyl chains are proportionally represented in both the mature cardiolipin profi le and the available pool of acyl chains. This is likely because the organs generally maintain their cardiolipin pattern after development, although cardiolipin remodeling continues during aging with only subtle changes to the acyl chain composition or distribution ( 16,(46)(47)(48).

Dynamic simulation of hepatic cardiolipin molecular species
Visualization of the distribution of predicted immature cardiolipin molecular species generated by the distribution of hepatic mitochondrial PG molecular species revealed a more diversifi ed profi le of species compared with those of myocardial immature cardiolipin species ( Fig. 4A ). Additionally, because molecular species from immature cardiolipin are minimally represented in the mature cardiolipin, only a minimal decay rate was used (between 0.00001 and 0.00005) for the dynamic simulation representing a fraction of de novo cardiolipin biosynthesis. The selection of these decay rates, however, do not address specifi c kinetic fl ux of the enzymes involved in synthesis and are only designed to refl ect the relative representation of de novo cardiolipin synthesis in the cardiolipin distribution. The acyl chain distribution of hepatic cardiolipin contained more linoleic acyl chains compared with that of myocardial cardiolipin. The seeded acyl chain distribution used for hepatic cardiolipin dynamic simulation consisted of: 16 ( Fig. 4C ). The simulation achieved an equilibrium around 250 iterations ( Fig. 4D ). This equilibrium state resembles the exact distribution for hepatic cardiolipin obtained from MDMS-SL analysis of the hepatic mitochondrial lipidome ( Fig. 4E ). The major and minor molecular species identifi ed by MDMS-SL also match the species generated by the dynamic simulation ( Table 2 and supplementary materials). Moreover, the level of detection of hepatic cardiolipin molecular species demonstrated by the dynamic simulation is similar to that of the heart (i.e., between 0.1% and 0.2%). Thus, the simulation strategically facilitates detection of low to very low abundance isomers of cardiolipin species, the levels of which are below the detection limitation of MDMS-SL analysis at its current stage. tetra-18:2 cardiolipin. In addition, we also observed several very low abundance fragments at m/z 253. 21 ( Fig. 6A, inset a ). This observation indicates the region corresponding to 18:2 carboxylate, along with a fragment at m/z 415.22 corresponding to the deprotonated 18:2 lysophosphatidic acid, and a product ion at m/z 592.36 corresponding to the loss of a 18:2 ketene ( Fig. 6A ). These fragment ions identifi ed the ion as a  ( Fig. 6B , inset d). Therefore, product ion analysis using an Orbitrap-LTQ mass spectrometer validated the accuracy of the dynamic simulation results that indicate the value of this approach for the predictive analysis of very low abundance cardiolipin isomeric species that may not be readily detectable by the current stage of MS. This identifi cation using product ion analysis is consistent with the results from the dynamic simulation ( Fig. 6A, inset b ).
By using this approach, we also validated the presence of very low abundance cardiolipin isomers as predicted by dynamic simulation residing at much lower abundance ions compared with the dominant ion at m/z 723. 48  Mouse hepatic mitochondrial lipids were extracted by a modifi ed Bligh and Dyer procedure as described in the section of "Materials and Methods." The cardiolipin molecular species in the lipid extracts were identifi ed and quantifi ed as previously described ( 15 ). The mass content and distribution of cardiolipin species represent the mean of four separately prepared mitochondria from different animals. The distribution of cardiolipin molecular species determined by the dynamic simulation represents the one frame output 0.1 ␣ 250 iterations for liver cardiolipin. This data was obtained using the MS output function. Values < 0.01 percent are not displayed but can be found in supplementary tables. CL, cardiolipin.
is commonly believed that acyl CoA as well as PC and PE contribute to the pool of acyl chains used for cardiolipin remodeling in mitochondria ( 20,21,23,24,38,43 ). In addition, Lands cycle remodeling utilizes the sn -2 positions of glycerophospholipids for remodeling ( 49,50 ).

Prediction of acyltransferase and transacylase activities as well as acyl selectivities in cardiolipin remodeling from different organs
Analysis of the mitochondrial lipidome from the heart, liver, and brain allowed for the opportunity to characterize possible substrates used for cardiolipin remodeling. It   Fig. 5. Dynamic simulation of brain cardiolipin molecular species. A: Predicted brain immature cardiolipin molecular species were generated representing at time point 0. Generated distributions of cardiolipin molecular species were represented at iteration 50 (B), 100 (C), and 250 (D). E: A representative mouse brain cardiolipin distribution as determined by MDMS-SL. The asterisks indicate the individual ion peaks of mouse brain cardiolipin. The mouse brain cardiolipin molecular species detected and quantifi ed by MDMS-SL as well as determined by the dynamic simulation were tabulated in Table 3 .  Mouse brain mitochondrial lipids were extracted by a modifi ed Bligh and Dyer procedure as described in the section of "Materials and Methods." The cardiolipin molecular species in the lipid extracts were identifi ed and quantifi ed as previously described ( 15 ). The mass content and distribution of cardiolipin species represent the mean of four separately prepared mitochondria from different animals. The distribution of cardiolipin molecular species determined by the dynamic simulation represents the one frame output 0.1 ␣ 250 iterations for brain cardiolipin. This data was obtained using the MS output function. Values < 0.01 percent are not displayed but can be found in supplementary tables. CL, cardiolipin. these enzymatic activities to the individual pools. It should be emphasized that there exist multiple acyltransferases and transacylases that have demonstrated cardiolipin remodeling activity, and indeed a variety of isoforms have been detected of these enzymes ( 20,38,(51)(52)(53). The current dynamic simulation approach is unable to differentiate these isoform activities yet proposes the relative selectivity of the combined action of various acyltransferases or transacylases. Additionally, we do not infer any selection toward remodeling specifi c de novo synthesized molecular species; instead, we infer that all immature cardiolipin species are remodeled equally and that partially and mature remodeled cardiolipin molecular species are treated with equal affi nity.
To rebuild the acyl chain distribution of cardiolipin, the overall distribution of acyl chains from the PC and PE sn -2 positions was calculated assuming no acyl selectivity originating from a transacylation reaction. This approach was followed determining the contribution of acyl CoA fatty acids with designated acyl selectivity for particular fatty acids. There are multiple reasons for taking this approach. First, we were attempting to determine if there was common acyl selectivity for acyl CoA between organs with dif-in molecular species at the sn -1 and sn -2 positions, allowing us to infer Lands cycle remodeling based on acyl chain identifi cation. Determination of the original sources (i.e., acyl CoA, PC sn -2, and PE sn -2 acyl chains, and PG molecular species) of the acyl pool of cardiolipin would facilitate the assessment of the enzymatic activities involved in cardiolipin biosynthesis and remodeling. In addition, modeling the molecular mechanism in which acyl chains are transferred into the sn -1, sn -1' or sn -2, sn -2' positions of cardiolipin either equally or with designated selection identifi es the potential components relegating the remodeling machinery. Furthermore, through rebuilding the acyl chain distribution that serves as the acyl pool in the remodeling process by proper selection of acyl chains from different sources, the dynamic simulation allows us to assess these activities. Specifi cally, the decay rate of the synthesized immature cardiolipin species calculates the de novo synthesis activity, the amounts of acyl chains contributed from PC and PE sn -2 fatty acyl pools indicate the transacylase activities, the amounts of acyl CoA used in the remodeling determine the acyltransferase activities, and the differential utilization of the different acyl chains from the PC, PE, or acyl CoA pool represents the selectivities of  Tables III and IV). ferent acyl CoA compositions as well as cardiolipin distributions. If similarities existed in acyl selectivity, then it would suggest a common enzyme or isoform of an enzyme likely involved in the remodeling process among tissues. Indeed, it is not unlikely that numerous acyltransferases exist and that the relative selectivity refl ects the coordinated activities of these acyltransferases; however, infer- Fig. 7. The contribution and selectivity of acyl chains for cardiolipin remodeling from different fatty acyl donor pools determined by dynamic simulation. Based on the acyl chains used to generate mature myocardial, hepatic, and brain cardiolipin, inference was used to determine the possible contribution from acyltransferase, transacylase, or cardiolipin biosynthesis activities. The possible proportion of acyl chains from PG, PC sn -2 positions, PE sn -2 positions, and acyl CoA as determined by MDMS-SL are represented in the chart that would match the acyl chain distribution used in the dynamic simulation. Selectivity of specifi c acyl chains for acyltransferase dependence (acyl CoA) (A) or transacylase dependence (PC or PE sn -2 acyl chains) (B) was inferred in the heart, liver, and brain. This inferences approach allows the comparison of different acyl selectivities for both acyltransferase and transacylase activities.
cardiolipin. In contrast, the proportional selectivity for 20:2 and 20:3 was only found in myocardial and hepatic cardiolipin and not brain cardiolipin. Additionally, there appeared to be selectivity for 16:1 acyl chains in acyl CoA for hepatic cardiolipin and to even a greater extent for brain cardiolipin; however, no selectivity for 16:1 could be concluded for myocardial cardiolipin. The similarities and differences in acyl selectivity between tissues may result from tissue-specifi c splicing variants or post-translation modifi cations of various known acyltransferases ( 20,38,41,45,52,54 ). Recently, the acyltransferase MLCL AT demonstrated a relative selectivity of 18:2 > 18:1 > 16:1 ( 40 ), which is consistent with the acyl CoA selectivity in the examined heart and liver tissue and also in the examined brain tissue except for 16:1 CoA selectivity.
We also modeled the acyl selectivity of transacylase activity ( Fig. 7B ). This approach demonstrated variability in the selectivity for 16:1, 18:1, 18:2, 20:2, and 20:3 acyl chains between different organs with no distinct pattern. Because it is likely that cardiolipin remodeling might be relegated by the coordinated activities of both acyltransferases and transacylases, the acyl selectivity from both activities might generate the acyl chain pool used for cardiolipin remodeling. In order to not infer any intrinsic bias in this analysis, we attempted to highlight the relative selectivities of acyl chains in solely acyltransferase or transacylase dependence in the current study. In all likelihood, the combination of the two mechanisms occur in synchrony to facilitate the highly specifi c acyl chain distributions present in different tissues. Additional experimentation might also reveal the relative proportion of acyltransferase versus transacylase activities contributing to the various cardiolipin distributions.

DISCUSSION
Recent advances in lipidomic analysis now allow for the accurate quantifi cation of hundreds to thousands of lipid molecular species. However, development of bioinformatics and algorithmic approaches to integrate lipidomic data into a systems biology interpretation will be required to transition lipidomics from a descriptive to a mechanistic paradigm. In this study, we utilized the MDMS-SL platform to structurally identify and quantify the lipid molecular species of heart, liver, brain, and lung mitochondria and then used this information to interpret the biological mechanisms governing cardiolipin remodeling. This involved the computational inferences of known biological parameters important in lipid remodeling. In addition, we developed a series of software tools that can be used to compare and visualize data obtained from MDMS-SL analysis as well as extend the range and detection of lipidomic analysis. This composite experimental and bioinformatic approach represents a major advance in interpreting the complex heterogeneity of the cellular lipidome.
Cardiolipin is intricately involved in maintaining bioenergetic effi ciency and homeostasis of cellular metabolism ( 11,13 ). The complex biophysical properties as well as the Alternatively, we attempted to fi t the acyl chain pool used for cardiolipin remodeling with a proportion of acyl chains from sn -2 positions of PC and PE as well as acyl CoA inferring some acyl chain selectivity. This approach allowed us to increase the percentage contribution of selected acyl chains as well as cardiolipin species from predicted immature cardiolipin if these species were visualized in the mature cardiolipin profi le. The distribution of acyl chains used for the myocardial cardiolipin simulation could be formed by contributing 4% and 3% each from PC and PE sn -2 positions, respectively, without acyl selectivity. The remaining 93% of the acyl chains would be donated from acyl CoA with the designated acyl specifi city ( Fig. 7A ). This specifi city was most evident in 18:2, 20:3, and 22:6 acyl chains and less selective for 18:1, 20:2, and 22:5 acyl chains.
Due to the lower amount of 20:4 and 22:6 containing acyl chains in hepatic cardiolipin, the level of transacylation from PC and PE appeared to be lower in the hepatic than in the myocardial cardiolipin. This was represented by 1% each for both glycerophospholipid classes without acyl selectivity. The remaining 98% of the acyl chains were contributed by acyl CoA with selectivity assigned to 18:2 and 20:3. Less selectivity was designated to 16:1, 18:1, 20:2, and 22:6 ( Fig. 7A ).
The composition of acyl chains in brain cardiolipin was more complex than that found in myocardial and hepatic cardiolipin. The most abundant acyl chains in brain cardiolipin are 18:1, 20:4, and 22:6, with lesser amounts of 18:2. These aforementioned acyl chains are abundant in the sn -2 position of PC and PE. This suggests that transacylation activity is greater for brain cardiolipin than for myocardial or hepatic cardiolipin. In addition, the proportion of saturated acyl chains was greater in brain cardiolipin than hepatic and myocardial cardiolipin, suggesting that immature cardiolipin contributes signifi cantly to the species composition of mature cardiolipin in brain than it does for hepatic and myocardial cardiolipin. Based on these criteria, it was determined that 10% of the acyl chains of brain cardiolipin molecular species likely originated from predicted immature cardiolipin, 22% and 14% of the acyl chains nonselectively originated from the sn -2 position of PC and PE, respectively, and the remaining 64% of the acyl chains originated from acyl CoA with acyl selectivity. The acyl chains that displayed the greatest selectivity were 18:2, 22:6, 16:1, and to a lesser extent 22:5 and 18:1 ( Fig. 7A ).
Analysis of the designated acyl selectivity in diverse acyl CoA distributions revealed an unexpected fi nding when comparing all three organs ( Fig. 7 ). Aligning the determined acyl selectivities for acyl CoA revealed the possibility that there are consistent selectivities among tissues as well as some differences. The relative proportional selectivity for 18:1 CoA and 18:2 CoA appeared to be consistent among all tissues. This was even more surprising for brain cardiolipin, which contains more 18:1 CoA and less 18:2 CoA content than that in other tissues. The proportional selectivity for acyl CoA in 22:5 and 22:6 acyl chains was evident in myocardial and brain cardiolipin, but not hepatic acyl chains used during the remodeling process. Furthermore, it provides precise predictions as to the relative composition of isomeric molecular species, which can clarify the identifi cation of low abundance molecular species of cardiolipin, which may not be detectable by current analytical technologies. This approach offers a systematic method of analyzing similarities between cardiolipin remodeling among different tissues, a key step toward dissecting the functions of the enzymes involved in lipid regulation. diverse molecular species of cardiolipin have stimulated the investigation of numerous analytical and mathematical approaches to determine the importance of cardiolipin in biological function. Coarse-grained molecular dynamics, thermodynamics and thermal response, atomic force microscopy, and atomic-scale simulations have all been used to characterize structural orientation and interaction of cardiolipin in membrane dynamics (55)(56)(57)(58)(59)(60). More recently, mathematical models using linear algebraic representation have been used to explain the acyl selectivity of cardiolipin remodeling based on free energies ( 61 ). Although that approach demonstrated novel insight into potential mechanisms governing acyl selectivity, it did not directly address the importance of kinetics in the generation of cardiolipin acyl chain distributions. In addition, that approach only addressed distributions that are relegated by tetralinoleic abundance of cardiolipin; however, brain cardiolipin does not conform to that model ( 16,26 ). Moreover, alterations in free energies are complex functions of enthalpy and entropy driven by altering molecular membrane symmetry and the order/disorder at the membrane interface.
A critical point addressed by the cardiolipin dynamic simulation in the current study was imposing the option for sn -1 or sn -2 selectivity in cardiolipin acyl chain remodeling. This option allowed for the inference of the rate as well as pattern distribution of cardiolipin molecular species if acyl chain selectivity was employed or if all four chains were remodeled equally to reach the equilibrium distribution. Our fi ndings indicate that treating all four chains equally for remodeling generated a cardiolipin distribution most representative of that obtained by MDMS-SL. This assumption supports the fi ndings of Malhotra et al. ( 62 ), who demonstrated that sn -1 and sn -2 positions of cardiolipin are remodeled equally using purifi ed recombinant Drosophila tafazzin.
The use of a dynamic simulation approach for identifying and modeling the biological parameters regulating cardiolipin synthesis and remodeling represents an important advance in integrating lipidomic data with the biochemical mechanisms of cardiolipin remodeling. This approach exploits differences in cardiolipin composition and metabolism between tissues ( 63 ) to gain mechanistic insight into cardiolipin remodeling pathways and kinetics. Similar approaches using algorithmic inference and in silico kinetic modeling to describe the enzymatic parameters governing the molecular species distribution found in other lipid classes are also possible through this approach. The integration of a fl ux-based approach for lipidomics has been utilized in describing eicosanoid metabolism; however, the number of parameters used to ascertain the fl ux of biological pathways can result in overfi tting the experimental data, depending on the model ( 64 ).
The dynamic simulation provides a formidable tool for analyzing the complex organization of cardiolipin molecular species as controlled by biosynthesis and remodeling. The strategy that we presented allows for the inference of acyltransferase and transacylase activities as well as acyl chain selectivities that would be needed to generate the