The genome-scale metabolic model (GEM) of C. acnes, iCA843, was reconstructed based on the established procedures (Thiele and Palsson, 2010). To reconstruct the GEM of a virulent C. acnes strain, the genome sequences and corresponding annotation of HL043PA1 from National Center for Biotechnology Information (NCBI) database as of 22 March 2022. HL043PA1 is a strain belonging to ribotype 5, which is known to be associated with acne vulgaris. The preliminary metabolic network draft was built using CarveMe (Machado et al., 2018). The nomenclature of metabolites and reactions was based on the BiGG database (Schellenberger et al., 2010). Next, additional metabolic and transport reactions of C. acnes were annotated by EggNOG 5.0 (Huerta-Cepas et al., 2019) and BlastKoala (Kanehisa et al., 2016), with various database such as BiGG (King et al., 2016), KEGG (Kanehisa, 2000), SEED (Henry et al., 2010), UniProt (Bateman et al., 2021), and TransportDB (Elbourne et al., 2017). Then, the relevant metabolic reactions were included with the corresponding gene-protein-reaction (GPR) assignments based on either direct or indirect biochemical evidence from the biochemical databases and literature. An effort was made to annotate metabolic gap-filled reactions that lacked GPR annotation, and in cases where GPR information was unavailable, the gap-filled reactions were excluded from the model. To enhance the reliability of each reaction in the model, a comparison of protein sequences using BLASTp (Altschul et al., 1990) was conducted prior to their inclusion with e-values< 1 × 10^-50 and a percentage identity >70%. The reversibility of coupled reactions was corrected based on the biochemical information from literature and online databases such as MetaCyc (Caspi et al., 2014), Brenda (Chang et al., 2021), Virtual Metabolic Human (VMH) (Noronha et al., 2019), and eQuilibrator (Beber et al., 2022), to comprehensively consider the physiological direction and biochemical thermodynamics of the reactions. Next, the model leveraged with previously reconstructed metabolic model of P. acnes 6609 (McCubbin et al., 2020) to cross-check the reactions with GPR associations and further expand the energy metabolism. To do this, homologous proteins between the HL043PA1 and 6609 strains were identified based on BLASTp, and the metabolic reactions associated with the resulting homologous proteins were inspected manually based on literature and gene annotation. Only reactions confirmed to have protein homology and literature evidence, and a valid GPR were added to the current reconstruction. For common reactions between the current and previous reconstruction, the GPR of these reactions was compared and updated accordingly if necessary. In particular, for differing reactions between the models, extensive research was conducted based on protein homologs, biochemical databases, and literature to determine whether to include these reactions, and they were included in the current model only if biological evidence was present. At the very last step, the quality of reconstructed GEM was evaluated by comparing it with relevant experimental data and a metabolic model test suite called MEMOTE (Lieven et al., 2020), which is a community standard for this purpose.

We derived a biomass equation for C. acnes with information of macromolecular and monomeric composition information obtained from published data (see Supplementary Data Sheet 1 ). Note that we also partially used biomass composition data from a taxonomically close species in the case where C. acnes specific data is unavailable (Rocha et al., 2008). Note that the macromolecular composition and some parts of the monomeric composition including lipid, polysaccharide and small molecules, were referred from the Propionibacterium biomass equation derived from the previous model (McCubbin et al., 2020), since taxonomically related species are shared among Propionibacterium. The protein, DNA and RNA compositions of C. acnes were estimated based on the genome sequence data used in the model reconstruction, while the fatty acid composition of C. acnes was obtained from literature (Moss et al., 1967). The growth and non-growth associated maintenance (GAM) were calculated based on the macromolecular composition, and non-growth associated maintenance (NGAM) were assumed to be identical to that of P. acidipropionici (Zhang et al., 2015).

The metabolic behavior of C. acnes was simulated under different conditions using constraint-based FBA (Raman and Chandra, 2009; Orth et al., 2010). All constraint-based flux analysis simulations were carried out using COBRA Toolbox in MATLAB R2020a (Schellenberger et al., 2011) with CPLEX optimization solver. The constraints used in model simulations are provided in Supplementary Table 1 . The cell was set to maximize biomass objective function ( v b i o m a s s ) while constraining uptake rates of other nutrients, such as uptake rates of carbon source and other complex medium components. The glucose and glycerol uptake rate were constrained at 10 mmol gDW^-1 h^-1 and 20 mmol gDW^-1 h^-1, respectively. Mathematical representation of the optimization problem can be expressed as follows:

Where S i j   is the stoichiometric coefficient of metabolite i that participates in reaction j and v j is the flux of reaction j. Reaction constraints were given by assigning lower and upper bounds to the reaction ( v j ) as v j m i n and v j m a x , respectively.

We predicted vitamin auxotrophy by constraining the uptake rate of each vitamin to zero and maximizing the biomass objective function. Vitamin auxotrophy was defined as over a 90% decrease in growth when the corresponding vitamin was excluded from the media (Koduru et al., 2017). SCFA production rate and intracellular flux distributions in glucose and glycerol were simulated using a variation of parsimonious FBA (pFBA) that optimizes the objective function and sequentially minimizes the flux through the model (Lewis et al., 2010). The turnover rates of metabolites under the aforementioned conditions were described using flux-sum, as in previous studies (Chung and Lee, 2009; Mishra et al., 2016). Assuming that the cell is in steady-state flux balanced condition, the flux-sum for metabolite i , ϕ i , is given by:

The sum of all fluxes containing metabolite i will give the total rate of consumption and generation, having this value will give the turnover rate of metabolite i . The flux and flux-sum intensity were also calculated in our study. Flux intensity was obtained by normalizing the flux values of each reaction by the maximum flux value of the reaction across the conditions, while flux-sum intensity for each reaction was obtained by dividing them by the corresponding maximum flux-sum value of the metabolite across the conditions. The flux and flux-sum values for both glucose and glycerol conditions are provided in Supplementary Table 2 .

As potential antimicrobial targets for C. acnes, essential genes of C. acnes in glycerol condition were identified by using the single gene deletion function provided in the COBRA toolbox (Joyce and Palsson, 2007). Genes whose knockouts resulted in a predicted growth rate of less than 10% of the wild-type predicted growth rate were considered essential. To sort out the essential genes in C. acnes that are homologous to other skin microbe gene sequences, the essential genes were further subjected to a BLASTp protein sequence similarity search against 180 abundant skin bacteria taxa present at > 0.1% of the reads in at least one sample among 251 collected samples (Bewick et al., 2019). The list of the abundant microbiome taxa is provided in Supplementary Table 3 . The whole genome sequences of reference strains of abundant skin microbiome taxa were retrieved from NCBI and used for BLASTp analysis. The protein sequences with an e-value < 1 × 10^-50 were considered homologs (Heinken et al., 2023).

Distinct associations have been observed between different ribotypes (RTs) of C. acnes and acne vulgaris (Fitz-Gibbon et al., 2013), indicating that exploring the metabolic differences among RTs may provide new insights into acne vulgaris. To investigate the variations in metabolic and phenotypic behavior among RTs, we focused on two specific RTs for GEM reconstruction and analysis: RT1, which is highly prevalent in both affected and healthy skin, and RT6, which is associated with healthy skin. The genome sequences of C. acnes strains ATCC6919 and HL110PA3 were utilized for the reconstruction of RT1 and RT6 GEMs, respectively. These genome sequences were retrieved from the NCBI database as of 23 February 2023.

The RT-specific GEMs were built based on protein homology and the reconstructed iCA843 model. SonicParanoid (Cosentino and Iwasaki, 2019) was utilized in its default mode to search for homologous proteins between ATCC6919 (RT1) and HL043PA1 (RT5), as well as between HL110PA3 (RT6) and HL043PA1 (RT5). The identified homologous proteins were used to transfer reactions with appropriate GPR associations from iCA843 to the RT1 or RT6 GEMs. Reactions lacking GPR information were automatically transferred to other GEMs if they were necessary for the growth or metabolic phenotypes of C. acnes. Reactions associated with RT-specific proteins were identified based on the draft GEMs reconstructed using ModelSEED (Seaver et al., 2021) and CarveMe (Machado et al., 2018) with manual curation.

To identify orthologous and ribotype-specific genes within the C. acnes strains, a total of 167,830 protein sequences was compiled from 73 strains representing different ribotypes. These sequences were then subjected to clustering using CD-HIT program (version 4.8.1) with an amino acid similarity threshold of 70% (Fu et al., 2012). From the pool of dispensable genes found in two or more strains, we specifically retained the genes that were unique to RT1, 5, or 6.

The GEM for C. acnes HL043PA1, iCA843, was reconstructed following procedures (see Materials and Methods). Initially, we build a draft metabolic network based on the functional annotation of genes from the whole genome sequence of the HL043PA1 strain. This strain belongs to RT5, which is reported to be strongly associated with acne vulgaris (Fitz-Gibbon et al., 2013; Barnard et al., 2016). The network was then manually checked to identify any discrepancies between the network and known physiological metabolism. In this step, it is important to rectify incorrect GPR relations by considering species-specific enzyme annotation and include the relevant metabolic reactions which may not be functionally assigned due to the limitation of homology-based approaches based on biochemical databases and literatures. As an example, we updated the information regarding methylmalonyl-CoA carboxyltransferases (MMC; EC 2.1.3.1), a critical reaction in the propionate biosynthesis. We made the necessary modification to the GPR annotation of MMC by newly incorporating a previously unannotated subunit that was identified as a hypothetical protein lacking specific functionality in the NCBI annotation. We also added type-2 phosphatidic acid phosphatase (EC 3.1.3.4), which was not included initially in the draft model, according to the literature reporting that it is required for de novo synthesis of phosphatidylglycerol in phospholipid metabolism (Kanoh et al., 1999). As a result, the GPR of 198 reactions were updated and 392 metabolic reactions were newly included in the model, based on the functional annotation, biochemical databases, and experimental evidence. In addition, a total of 705 metabolic reactions were excluded from the draft model due to redundancy in fatty acid and phospholipid metabolism, as well as the presence of periplasmic reactions that are not feasible in Gram-positive bacteria such as C. acnes. Next, the reconstructed model was further expanded by integrating existing metabolic model of P. acnes 6609 (McCubbin et al., 2020), which was built based on the pan-genome of Propionibacterium and offered novel metabolic insights into the energy conservation mechanism. This process led to an update of GPR for 297 existing reactions, deletion of 94 misannotated reactions, and incorporation of 163 new reactions, including a novel ferredoxin-based energy conservation mechanism of Propionibacterium proposed in previous study. Reversibility of all reactions were cross-checked based on several databases such as MetaCyc (Caspi et al., 2014), Virtual Metabolic Human (VMH) (Noronha et al., 2019), and eQuilibrator (Beber et al., 2022), which prevents biologically unfavorable intracellular fluxes in model simulation. In addition, we formulated a biomass equation which was derived from the macromolecular and monomer compositions as described in Material and Methods, and Supplementary Data Sheet 1 . Finally, the resulting model, iCA843, comprises 843 genes, 1510 reactions and 1194 metabolites, covering comprehensive central carbon, amino acid, and lipid metabolisms, as well as relevant biosynthetic pathways of key virulence factors such as SCFAs, TAG lipase and various types of porphyrins ( Figure 1 ). The list of reactions, along with their respective gene associations and metabolites, is provided in both systems biology markup language (SBML) and excel format ( Supplementary Data Sheet 2 ). Additionally, we confirmed the network consistency of iCA843 using the online tool MEMOTE (Lieven et al., 2020), achieving an overall score of 89% ( Supplementary Data Sheet 3 ).

Next, we assessed the quality of iCA843 by simulating the cell growth under glucose condition. Remarkably, the predicted fermentative behavior was highly consistent with the reported SCFA measurements (Stowers et al., 2014), with propionate producing more than twice that of acetate ( Figure 2A ). It should be noted that the secretion of propionate and acetate was not constrained in the FBA simulation (see Materials and Methods), showing that iCA843 successfully captures the physiological and metabolic traits of C. acnes, a representative propionic acid bacteria (PAB) that efficiently ferments carbon sources to produce propionate through the Wood-Werkman cycle. The Wood-Werkman cycle consists of MMC, propionyl-CoA:succinate CoA transferase, and methylmalony-CoA mutase (Bücher et al., 2021). Next, we performed simulations to determine the essential vitamins B and C required for cell growth, which identified 4 auxotrophic vitamins (thiamin, riboflavin, pantothenate, and cobalamin) and 5 prototrophic vitamins (nicotinamide, pyridoxine, biotin, folate, and ascorbate) ( Figure 2B ). Notably, our results are consistent with experimental observations for 3 vitamins, namely nicotinamide, pantothenate, and biotin (McDowell et al., 2016), while there is a discrepancy between our simulation results (thiamin as essential for growth) and the experimental evidence. This inconsistency suggests the necessity for further investigation to reconcile the differences and gain a deeper understanding of the role of thiamin in the growth of C. acnes. Upon conducting a protein sequence similarity search targeting known enzyme sequences for thiamin biosynthesis, a complete biosynthesis pathway for thiamin could not be identified, making it impossible to gap-fill the related pathway. This finding suggests the presence of potential knowledge gaps in the genome annotation or an alternative, undiscovered biosynthesis pathway, which may explain the HL043PA1-specific auxotrophic behavior. In addition, we conducted fermentable substrate phenotyping for various carbon sources, and compared in silico predictions with experimental data reported by Puhvel (1968), giving rise to the growth phenotypes which are in good agreement with 13 out of 14 different carbon sources available naturally under human skin environment or provided from skin care products ( Figure 2C ).

To provide an overview of C. acnes metabolic traits as a skin inhabitant and to highlight the uniqueness of C. acnes metabolism in relation to acne vulgaris, we compared iCA843 with GEMs of other skin bacteria that have been reported to contribute to acne vulgaris pathogenesis, Staphylococcus epidermidis and Klebsiella pneumoniae ( Figure 2D ) (Kumar et al., 2016; Henry et al., 2017; Díaz Calvo et al., 2022). All three species can utilize glycerol, a carbon source abundantly available in the skin environment that provokes bacterial fermentation (Balasubramaniam et al., 2020). They also commonly produce protoporphyrin IX, which is a precursor to heme, as well as several SCFAs such as acetate and butyrate. However, in contrast, C. acnes can uniquely synthesize coproporphyrin III, which is more relevant to acne lesions than protoporphyrin IX (Patwardhan et al., 2017), and utilize distinctive SCFA biosynthetic pathways to produce propionate and acetate as major metabolic byproducts through the Wood-Werkman cycle and the newly recognized succinyl-CoA:acetate CoA-transferase (SUCOAACTr) (Zhang et al., 2021), respectively. On the other hand, K. pneumoniae produce propionate from the propanediol pathway (Luo et al., 2012), while S. epidermidis has no metabolism for its synthesis. Furthermore, an analysis of the genome and model characteristics of the three skin pathogens revealed that although the C. acnes model has the shortest genome length and the fewest open reading frames, it also has the highest percentage of accounted ORFs in the model and gene-associated reactions ( Figure 2E ). In addition, a comparison between the genome and model characteristics of iCA843 and the previous P. acnes 6609 GEM showed that iCA843 encompasses a greater proportion of open reading frames and gene-associated reactions. The inclusion of a higher percentage of gene-supported reactions and a more comprehensive set of network components in iCA843 has significantly expanded the metabolic capabilities of the model. This expansion is expected to capture specific metabolic characteristics unique to C. acnes.

C. acnes exhibits a significant fermentative characteristic by producing a substantial quantity of propionate through the Wood-Werkman cycle. This ATP-independent pathway efficiently converts pyruvate into oxaloacetate (Wang et al., 2015). In addition to propionate, C. acnes also produces other SCFAs, including acetate. These SCFAs play a role in stimulating free fatty acid receptors and can trigger inflammatory reactions in skin immune cells (Sanford et al., 2019), while propionate can have cytotoxic effects by causing pH changes (Tax et al., 2016), and can also provoke an immune response (Kim et al., 2002; Nagy et al., 2005). The production of SCFAs by C. acnes contributes to the complex interplay between the microbiota and the host immune system in the context of skin health and inflammation. Therefore, our aim is to elucidate the metabolic flux distributions attaining to these fermentative behaviors in human skin, as an effort to gain a better understanding of the underlying intercellular mechanisms involved in the development of acne vulgaris. We performed flux simulations under anaerobic conditions to mimic the environment of the obstructed pilosebaceous unit in acne vulgaris (Shannon, 2020). Furthermore, we compared a comparison of the growth and metabolic state of C. acnes under two different conditions: glucose and glycerol. Glucose was considered as the control condition, while glycerol represented the major carbon source in the pilosebaceous unit (see Materials and Methods). This is due to the fact that human sebum primarily consists of triglycerides, fatty acids, squalene, and wax esters (Akaza et al., 2014), and C. acnes produces extracellular lipase enzymes that hydrolyzes the triglycerides present in sebum, leading to the release of glycerol as a nutrient source (Coenye et al., 2021). However, there is a lack of experimental evidence regarding the uptake of other sebum constituents by C. acnes. Therefore, glycerol is considered the primary endogenous carbon sources in human sebum, which aligns with the existing knowledge that glycerol serves as a major carbon source for the skin microbiome, facilitating growth and biosurfactant production (Saikia et al., 2012; Balasubramaniam et al., 2020). The simulation results showed that the growth rate under glycerol condition was 31.5% lower compared to glucose condition, and the production of acetate was negligible, while there was a 1.6-folds increase in propionate production rates observed in the glycerol condition ( Figure 3A ). These findings consistent with observations made in several Propionibacterium species (Liu et al., 2011; Wang and Yang, 2013; Zhang et al., 2015). The breakdown of carbon output from each source revealed that propionate exhibited the highest efflux in both conditions, while the contribution of CO₂ to the total carbon output was lower in the glucose condition as expected ( Figure 3B ).

The resulting internal fluxes within the central metabolism and SCFA biosynthetic pathways clearly revealed that pyruvate kinase (PYK) in the last step of glycolysis and all reactions in the Wood-Werkman cycle have significantly higher fluxes in the glycerol condition compared to the glucose (up to 2.5-folds), while pyruvate dehydrogenase (PDH) is highly active in the glucose condition (up to 7.2-folds), resulting in higher acetate secretion coupled with CO₂ production ( Figure 3C ). We further explored the Wood-Werkman cycle and PDH from a redox balance perspective by quantifying the turnover rates of energy cofactors such as ATP and NAD based on their flux-sum intensity values (see Materials and Methods). It should be noted that flux-sum can represent the metabolite pool size by summing up all incoming or outgoing fluxes associated with the metabolite (Chung and Lee, 2009). As a result, we observed a higher turnover rate of NAD (42.2%) in the glycerol-rich condition ( Figure 3D ), which is mainly contributed by glycerol-3-phosphate dehydrogenase (G3PD1) and malate dehydrogenase (MDH) reactions in glycolysis and TCA cycle, respectively. This observation is in good agreement with previous study on several Propionibacterium species, for example P. jensenii, (Liu et al., 2011; Wang and Yang, 2013; Zhang et al., 2015), which also possess the Wood-Werkman cycle to maintain cellular redox balance (Luna-Flores et al., 2018; Gonzalez-Garcia et al., 2020). Similarly, ATP turnover rate was increased by 20.4% in the glycerol due to the high ATP demand via glycerol kinase (GLYK) to consume the carbon source. Consequently, the cellular metabolism shifted towards ATP regeneration rather than cell growth, as indicated by the activation of PYK in glycolysis and the subsequent increase in the pyruvate pool which serves as a precursor for the metabolic reactions of the Wood-Werkman cycle. Overall, the simulation results suggest that in the context of acne vulgaris-associated skin conditions, C. acnes may utilize the Wood-Werkman cycle to replenish depleted NAD. This metabolic adaptation can lead to the overproduction of propionate, which in turn may trigger an inflammatory response in human skin.

Considering that the current antibiotics used for acne vulgaris treatment that reduce C. acnes population may cause dysbiosis in the skin microbiome, which can lead to other skin diseases (Chien et al., 2019; Thompson et al., 2020), we presented a model-driven framework for systematically screening antimicrobial candidates and identifying promising targets which selectively suppress the growth of C. acnes with minimal effects on other skin microbiota ( Figure 4A ). Initially, using iCA843, we applied gene essentiality analysis to determine the genes which are crucial for the cell growth (see Materials and Methods), and found 117 essential genes (13.9%) out of 843 genes (step 1). With the list of abundant microbiome taxa (180 species) (Bewick et al., 2019) and their whole genome sequences collected from NCBI database, the number of species containing homologous genes given each essential gene for C. acnes were obtained via protein sequence similarity search using BLASTp (step 2). It is followed by narrowing down the gene list which are found only within less than 5% and 1% of abundant skin microbes in step 3, resulting in 23 ‘unique’ and 3 ‘highly unique’ candidates, respectively ( Supplementary Table 5 ). In Figure 4B , we showed the distribution of essential genes including C. acnes specific antimicrobial targets (117 genes in total) in various metabolic subsystems which are classified based on their cluster of orthologous groups (COG) functional categories. The largest portion belongs to ‘coenzyme metabolism’ with two highly unique genes, indicating its high rigidity. Amino acids and nucleotide metabolisms have high number of essential genes, but most of them were not unique as expected since they are functionally conserved among bacterial species (Peregrín-Alvarez et al., 2009). Interestingly, 50% of the essential genes in the ‘lipid metabolism’ were found to be unique genes.

We further investigated biological mechanisms and pathways related to highly unique genes, which allowed us to identify promising C. acnes drug targets for the acne vulgaris treatment ( Figure 4C ). It should be noted that none of the highly unique candidates were homolog to human genome, and one of three genes have been discarded due to their functional ambiguity. The two highly unique genes encode dihydroneopterin aldolase (DHNPA2r) and 2-amino-4-hydroxy-6-hydroxymethyldihydropteridine diphosphokinase (HPPK2) which are involved in biosynthesis of tetrahydrofolate, a central cofactor in bacterial amino acid and nucleic acid metabolism (Tjong et al., 2022). In fact, a formate-tetrahydrofolate ligase enzyme, utilizing tetrahydrofolate as a substrate and producing 10-formyltetrahydrofolate, is present as a housekeeping gene in 72 strains of C. acnes (Kilian et al., 2012), indicating the essentiality of tetrahydrofolate for their survival. Thus, the current model-guided framework enabled us to identify two promising C. acnes-specific antimicrobial targets, the enzyme encoding DHNPA2r and HPPK2, which await further experimental validation.