Molecular determinants of protein evolvability.

The plethora of biological functions that sustain life is rooted in the remarkable evolvability of proteins. An emerging view highlights the importance of a protein's initial state in dictating evolutionary success. A deeper comprehension of the mechanisms that govern the evolvability of these initial states can provide invaluable insights into protein evolution. In this review, we describe several molecular determinants of protein evolvability, unveiled by experimental evolution and ancestral sequence reconstruction studies. We further discuss how genetic variation and epistasis can promote or constrain functional innovation and suggest putative underlying mechanisms. By establishing a clear framework for these determinants, we provide potential indicators enabling the forecast of suitable evolutionary starting points and delineate molecular mechanisms in need of deeper exploration.

Integrative single-cell transcriptomic analyses reveal the cellular ontological and functional heterogeneities of primary and metastatic liver tumors.

BACKGROUND: The global cellular landscape of the tumor microenvironment (TME) combining primary and metastatic liver tumors has not been comprehensively characterized. METHODS: Based on the scRNA-seq and spatial transcriptomic data of non-tumor liver tissues (NTs), primary liver tumors (PTs) and metastatic liver tumors (MTs), we performed the tissue preference, trajectory reconstruction, transcription factor activity inference, cell-cell interaction and cellular deconvolution analyses to construct a comprehensive cellular landscape of liver tumors. RESULTS: Our analyses depicted the heterogeneous cellular ecosystems in NTs, PTs and MTs. The activated memory B cells and effector T cells were shown to gradually shift to inhibitory B cells, regulatory or exhausted T cells in liver tumors, especially in MTs. Among them, we characterized a unique group of TCF7+ CD8+ memory T cells specifically enriched in MTs that could differentiate into exhausted T cells likely driven by the p38 MAPK signaling. With regard to myeloid cells, the liver-resident macrophages and inflammatory monocyte/macrophages were markedly replaced by tumor-associated macrophages (TAMs), with TREM2+ and UBE2C+ TAMs enriched in PTs, while SPP1+ and WDR45B+ TAMs in MTs. We further showed that the newly identified WDR45B+ TAMs exhibit an M2-like polarization and are associated with adverse prognosis in patients with liver metastases. Additionally, we addressed that endothelial cells display higher immune tolerance and angiogenesis capacity, and provided evidence for the source of the mesenchymal transformation of fibroblasts in tumors. Finally, the malignant hepatocytes and fibroblasts were prioritized as the pivotal cell populations in shaping the microenvironments of PTs and MTs, respectively. Notably, validation analyses by using spatial or bulk transcriptomic data in clinical cohorts concordantly emphasized the clinical significance of these findings. CONCLUSIONS: This study defines the ontological and functional heterogeneities in cellular ecosystems of primary and metastatic liver tumors, providing a foundation for future investigation of the underlying cellular mechanisms.

GPCRs and fibroblast heterogeneity in fibroblast-associated diseases.

G protein-coupled receptors (GPCRs) are the largest and most diverse class of signaling receptors. GPCRs regulate many functions in the human body and have earned the title of "most targeted receptors". About one-third of the commercially available drugs for various diseases target the GPCRs. Fibroblasts lay the architectural skeleton of the body, and play a key role in supporting the growth, maintenance, and repair of almost all tissues by responding to the cellular cues via diverse and intricate GPCR signaling pathways. This review discusses the dynamic architecture of the GPCRs and their intertwined signaling in pathological conditions such as idiopathic pulmonary fibrosis, cardiac fibrosis, pancreatic fibrosis, hepatic fibrosis, and cancer as opposed to the GPCR signaling of fibroblasts in physiological conditions. Understanding the dynamics of GPCR signaling in fibroblasts with disease progression can help in the recognition of the complex interplay of different GPCR subtypes in fibroblast-mediated diseases. This review highlights the importance of designing and adaptation of next-generation strategies such as GPCR-omics, focused target identification, polypharmacology, and effective personalized medicine approaches to achieve better therapeutic outcomes for fibrosis and fibrosis associated malignancies.

Apoptotic cell signals and heterogeneity in macrophage function: Fine-tuning for a healthy liver.

Functional heterogeneity in tissue macrophage populations has often been traced to developmental and spatial cues. Upon tissue damage, macrophages are exposed to soluble mediators secreted by activated cells, which shape their polarisation. Interestingly, macrophages are concomitantly exposed to a variety of different dying cells, which carry miscellaneous signals and that need to be recognised and promptly up-taken by professional phagocytes. This review discusses how differences in the nature of the dying cells, like their morphological and biochemical features as well as the specificity of phagocytic receptor usage on macrophages, might contribute to the transcriptional and functional heterogeneity observed in phagocytic cells in the tissue.

High-Complexity cellular barcoding and clonal tracing reveals stochastic and deterministic parameters of radiation resistance.

It is elusive whether clonal selection of tumor cells in response to ionizing radiation (IR) is a deterministic or stochastic process. With high resolution clonal barcoding and tracking of over 400 000 HNSCC patient-derived tumor cells the clonal dynamics of tumor cells in response to IR was analyzed. Fractionated IR induced a strong selective pressure for clonal reduction which significantly exceeded uniform clonal survival probabilities indicative for a strong clone-to-clone difference within tumor cell lines. IR induced clonal reduction affected the majority of tumor cells ranging between 96% and 75% and correlated to the degree of radiation sensitivity. Survival to IR is driven by a deterministic clonal selection of a smaller population which commonly survives radiation, while increased clonogenic capacity is a result of clonal competition of cells which have been selected stochastically. A 2-fold increase in radiation resistance results in a 4-fold (P < .05) higher deterministic clonal selection showing that the ratio of these parameters is amenable to radiation sensitivity which correlates to prognostic biomarkers of HNSCC. Evidence for the existence of a rare subpopulation with an intrinsically radiation resistant phenotype commonly surviving IR was found at a frequency of 0.6% to 3.3% (P < .001, FDR 3%). With cellular barcoding we introduce a novel functional heterogeneity associated qualitative readout for tracking dynamics of clonogenic survival in response to radiation. This enables the quantification of intrinsically radiation resistant tumor cells from patient samples and reveals the contribution of stochastic and deterministic clonal selection processes in response to IR.

Understanding natural killer cell biology from a single cell perspective.

During the last decade, advances in single cell technologies have ignited increased understanding of natural killer cells (NK cells), which turned out to be far more complex than originally thought. Ample studies have established tissue-specific phenotypic variation within this cell population; however, the functional implication of this vast variation is still unclear. At single-cell level, the function of a NK cell is tightly regulated by several checkpoints however upon proper recognition the cell can deliver a lytic hit as early as 10 min or could take hours before they can kill their target cells. Moreover, only a fraction of NK cells appears to kill target cells while the larger portion of NK cells appear to be non-cytotoxic. All these studies showed that the NK cell compartment is composed of cells with different functional strengths and efficacies, thereby highlighting the necessity of analytical platforms that allow the study of these important innate immune cells at single-cell level. In this review, we discuss and provide an overview on phenotypical and functional heterogeneity within the NK cell population and subsequently provide information regarding emerging technologies that highlight the importance of single-cell studies to understand the biology of these cells.

Functional Heterogeneity of Bone Marrow Mesenchymal Stem Cell Subpopulations in Physiology and Pathology.

Bone marrow mesenchymal stem cells (BMSCs) are multi-potent cell populations and are capable of maintaining bone and body homeostasis. The stemness and potential therapeutic effect of BMSCs have been explored extensively in recent years. However, diverse cell surface antigens and complex gene expression of BMSCs have indicated that BMSCs represent heterogeneous populations, and the natural characteristics of BMSCs make it difficult to identify the specific subpopulations in pathological processes which are often obscured by bulk analysis of the total BMSCs. Meanwhile, the therapeutic effect of total BMSCs is often less effective partly due to their heterogeneity. Therefore, understanding the functional heterogeneity of the BMSC subpopulations under different physiological and pathological conditions could have major ramifications for global health. Here, we summarize the recent progress of functional heterogeneity of BMSC subpopulations in physiology and pathology. Targeting tissue-resident single BMSC subpopulation offers a potentially innovative therapeutic strategy and improves BMSC effectiveness in clinical application.

Regionally encoded functional heterogeneity of astrocytes in health and disease: A perspective.

Increasing evidence has suggested that astrocytes demonstrate striking regionally allocated functional heterogeneity. Here, we discuss how this spatiotemporally encoded diversity determines the astrocytic phenotype along a finely grained spectrum from neuroprotective to deleterious states. With increasing recognition of their diverse and evolving roles in the central neuraxis, astrocytes now represent a tractable cellular target for therapies aiming to restore neural circuit integrity in a broad range of neurodegenerative disorders. Understanding the determinants of astrocyte physiology along with the true extent of heterogeneity in their regional and subregional functions will ultimately inform therapeutic strategy in neurodegenerative diseases.

Multiplicity and Polymorphism of Fish Hemoglobins.

The diversity of fish hemoglobins and the association with oxygen availability and physiological requirements during the life cycle has attracted scientists since the first report on multiple hemoglobin in fishes (Buhler and Shanks 1959). The functional heterogeneity of the fish hemoglobins enables many species to tolerate hypoxic conditions and exhausting swimming, but also to maintain the gas pressure in the swim bladder at large depths. The hemoglobin repertoire has further increased in various species displaying polymorphic hemoglobin variants differing in oxygen binding properties. The multiplicity of fish hemoglobins as particularly found in the tetraploid salmonids strongly contrasts with the complete loss of hemoglobins in Antarctic icefishes and illustrates the adaptive radiation in the oxygen transport of this successful vertebrate group.

The CD4-CD8- MAIT cell subpopulation is a functionally distinct subset developmentally related to the main CD8+ MAIT cell pool.

Mucosa-associated invariant T (MAIT) cells are unconventional innate-like T cells that recognize microbial riboflavin metabolites presented by the MHC class I-like protein MR1. Human MAIT cells predominantly express the CD8α coreceptor (CD8+), with a smaller subset lacking both CD4 and CD8 (double-negative, DN). However, it is unclear if these two MAIT cell subpopulations distinguished by CD8α represent functionally distinct subsets. Here, we show that the two MAIT cell subsets express divergent transcriptional programs and distinct patterns of classic T cell transcription factors. Furthermore, CD8+ MAIT cells have higher levels of receptors for IL-12 and IL-18, as well as of the activating receptors CD2, CD9, and NKG2D, and display superior functionality following stimulation with riboflavin-autotrophic as well as riboflavin-auxotrophic bacterial strains. DN MAIT cells display higher RORγt/T-bet ratio, and express less IFN-γ and more IL-17. Furthermore, the DN subset displays enrichment of an apoptosis gene signature and higher propensity for activation-induced apoptosis. During development in human fetal tissues, DN MAIT cells are more mature and accumulate over gestational time with reciprocal contraction of the CD8+ subset. Analysis of the T cell receptor repertoire reveals higher diversity in CD8+ MAIT cells than in DN MAIT cells. Finally, chronic T cell receptor stimulation of CD8+ MAIT cells in an in vitro culture system supports the accumulation and maintenance of the DN subpopulation. These findings define human CD8+ and DN MAIT cells as functionally distinct subsets and indicate a derivative developmental relationship.

Functional properties of single motor units in the human medial pterygoid muscle: Thresholds.

BACKGROUND: Little is known regarding the functional properties of single motor units (SMUs) in the medial pterygoid muscle (MPt) during jaw movements. OBJECTIVES: The aims are (a) to report the thresholds of onset of MPt SMUs during 4 goal-directed jaw movement tasks, and (b) to determine whether the threshold of onset of SMU activation varies with the velocity of jaw movement and the location within the muscle. METHODS: Intra-muscular electrodes were inserted in the right MPt of 18 participants performing ipsilateral (right), contralateral, protrusive and opening-closing jaw movements recorded at 2 velocities. Task phases were as follows: BEFORE, OUT, HOLDING, RETURN and AFTER. SMU onset thresholds were determined from the displacement (mm) of the lower mid-incisor point. Electrode location within 4 arbitrary muscle divisions was determined with computer tomography. Statistical tests: Spearman's correlations, Kruskal-Wallis tests; significance accepted at P < .05. RESULTS: A significant inverse relation occurred between velocity and threshold for the RETURN of the ipsilateral movement (n = 62 SMU thresholds), while a significant positive relation occurred for the OUT of the contralateral movement (n = 208); there were no significant associations for the protrusive (n = 131) and opening-closing (n = 58) tasks. Significant threshold differences occurred across the 4 muscle divisions only during the OUT of the contralateral and protrusive movements. Some evidence was provided for gender differences in MPt SMU properties. CONCLUSIONS: The absence of a significant inverse relation between velocity and SMU threshold for most recorded movements suggests the MPt acts as a stabilizer of the jaw in horizontal and opening-closing jaw movements.

Spatially Enriched Paralog Rearrangements Argue Functionally Diverse Ribosomes Arise during Cold Acclimation in Arabidopsis.

Ribosome biogenesis is essential for plants to successfully acclimate to low temperature. Without dedicated steps supervising the 60S large subunits (LSUs) maturation in the cytosol, e.g., Rei-like (REIL) factors, plants fail to accumulate dry weight and fail to grow at suboptimal low temperatures. Around REIL, the final 60S cytosolic maturation steps include proofreading and assembly of functional ribosomal centers such as the polypeptide exit tunnel and the P-Stalk, respectively. In consequence, these ribosomal substructures and their assembly, especially during low temperatures, might be changed and provoke the need for dedicated quality controls. To test this, we blocked ribosome maturation during cold acclimation using two independent reil double mutant genotypes and tested changes in their ribosomal proteomes. Additionally, we normalized our mutant datasets using as a blank the cold responsiveness of a wild-type Arabidopsis genotype. This allowed us to neglect any reil-specific effects that may happen due to the presence or absence of the factor during LSU cytosolic maturation, thus allowing us to test for cold-induced changes that happen in the early nucleolar biogenesis. As a result, we report that cold acclimation triggers a reprogramming in the structural ribosomal proteome. The reprogramming alters the abundance of specific RP families and/or paralogs in non-translational LSU and translational polysome fractions, a phenomenon known as substoichiometry. Next, we tested whether the cold-substoichiometry was spatially confined to specific regions of the complex. In terms of RP proteoforms, we report that remodeling of ribosomes after a cold stimulus is significantly constrained to the polypeptide exit tunnel (PET), i.e., REIL factor binding and functional site. In terms of RP transcripts, cold acclimation induces changes in RP families or paralogs that are significantly constrained to the P-Stalk and the ribosomal head. The three modulated substructures represent possible targets of mechanisms that may constrain translation by controlled ribosome heterogeneity. We propose that non-random ribosome heterogeneity controlled by specialized biogenesis mechanisms may contribute to a preferential or ultimately even rigorous selection of transcripts needed for rapid proteome shifts and successful acclimation.

Interprofessional teamwork in primary care: the effect of functional heterogeneity on performance and the role of leadership.

This study aimed to unravel the complexity of interprofessional teamwork in primary care teams by testing the relationship between functional heterogeneity and team performance through the mediating role of information elaboration, and the moderating roles of directive leadership and participative leadership. The moderated mediation model was validated using survey data from 1105 professionals and 97 supervisors in 143 Dutch primary care teams. The results confirmed the model and showed a significant negative effect of functional heterogeneity on information elaboration, which in turn had a positive effect on team performance. Both directive and participative leadership moderated the negative effect of functional heterogeneity on information elaboration to the extent that the indirect negative effect of functional heterogeneity on team performance became insignificant under high levels of either directive or participative leadership. The theoretical implications of these findings for the literature on healthcare, team diversity, and leadership, as well as the practical implications for policy makers, educationalists and managers of primary care teams, are discussed.

Single-cell analyses reveal functional classification of dendritic cells and their potential roles in inflammatory disease.

Dendritic cells (DCs) have crucial roles in immune-related diseases. However, it is difficult to explore DCs because of their rareness and heterogeneity. Although previous studies had been performed to detect the phenotypic characteristics of DC populations, the functional diversity has been ignored. Using a combination of flow cytometry, single-cell quantitative PCR, and bioinformatic analysis, we depicted the DC panorama with not only phenotypic but also functional markers. Functional classification of DCs in mouse lymphoid tissue (spleen) and nonlymphoid tissue (liver) was performed. The results revealed that expression of macrophage scavenger receptor 1 ( MSR1) and C-C motif chemokine receptors ( CCR) 1, CCR2, and CCR4 were elevated in liver DCs, suggesting increased lipid uptake and migration abilities. The enriched expression of costimulatory molecule CD80, TLR9, and TLR adaptor MYD88 in spleen DCs indicated a more-mature phenotype, enhanced pathogen recognition, and T-cell stimulation abilities. Furthermore, we compared DCs in the atherosclerotic mouse models with healthy controls. In addition to the quantitative increase in DCs in the liver and spleen of the apolipoprotein E-knockout ( ApoE-/-) mice, the functional expression patterns of the DCs also changed at the single-cell level. These results promote our understanding of the participation of DCs in inflammatory diseases and have potential applications in DC clinical assessment.-Shi, Q., Zhuang, F., Liu, J.-T., Li, N., Chen, Y.-X., Su, X.-B., Yao, A.-H., Yao, Q.-P., Han, Y., Li, S.-S., Qi, Y.-X., Jiang, Z.-L. Single-cell analyses reveal functional classification of dendritic cells and their potential roles in inflammatory disease.

Human papillomavirus oncoproteins induce a reorganization of epithelial-associated γδ T cells promoting tumor formation.

It has been shown that γδ T cells protect against the formation of squamous cell carcinoma (SCC) in several models. However, the role of γδ T cells in human papillomavirus (HPV)-associated uterine cervical SCC, the third-leading cause of death by cancer in women, is unknown. Here, we investigated the impact of γδ T cells in a transgenic mouse model of carcinogenesis induced by HPV16 oncoproteins. Surprisingly, γδ T cells promoted the development of HPV16 oncoprotein-induced lesions. HPV16 oncoproteins induced a decrease in epidermal Skint1 expression and the associated antitumor Vγ5+ γδ T cells, which were replaced by γδ T-cell subsets (mainly Vγ6+ γδlowCCR2+CCR6-) actively producing IL-17A. Consistent with a proangiogenic role, γδ T cells promoted the formation of blood vessels in the dermis underlying the HPV-induced lesions. In human cervical biopsies, IL-17A+ γδ T cells could only be observed at the cancer stage (SCC), where HPV oncoproteins are highly expressed, supporting the clinical relevance of our observations in mice. Overall, our results suggest that HPV16 oncoproteins induce a reorganization of the local epithelial-associated γδ T-cell subpopulations, thereby promoting angiogenesis and cancer development.

Functional Heterogeneity of Protein Kinase A Activation in Multipotent Stromal Cells.

Multipotent stromal cells (MSC) demonstrate remarkable functional heterogeneity; however, its molecular mechanisms remain largely obscure. In this study, we explored MSC response to hormones, which activate Gs-protein / cyclic AMP (cAMP) / protein kinase A (PKA) dependent signaling, at the single cell level using genetically encoded biosensor PKA-Spark. For the first time, we demonstrated that about half of cultured MSCs are not able to activate the cAMP/PKA pathway, possibly due to the limited availability of adenylyl cyclases. Using this approach, we showed that MSC subpopulations responding to various hormones largely overlapped, and the share of responding cells did not exceed 40%. Using clonal analysis, we showed that signaling heterogeneity of MSC could be formed de novo within 2 weeks.

Heterogeneous yet stable Vδ2(+) T-cell profiles define distinct cytotoxic effector potentials in healthy human individuals.

Human γδ T cells display potent responses to pathogens and malignancies. Of particular interest are those expressing a γδ T-cell receptor (TCR) incorporating TCRδ-chain variable-region-2 [Vδ2(+)], which are activated by pathogen-derived phosphoantigens (pAgs), or host-derived pAgs that accumulate in transformed cells or in cells exposed to aminobisphosphonates. Once activated, Vδ2(+) T cells exhibit multiple effector functions that have made them attractive candidates for immunotherapy. Despite this, clinical trials have reported mixed patient responses, highlighting a need for better understanding of Vδ2(+) T-cell biology. Here, we reveal previously unappreciated functional heterogeneity between the Vδ2(+) T-cell compartments of 63 healthy individuals. In this cohort, we identify distinct "Vδ2 profiles" that are stable over time; that do not correlate with age, gender, or history of phosphoantigen activation; and that develop after leaving the thymus. Multiple analyses suggest these Vδ2 profiles consist of variable proportions of two dominant but contrasting Vδ2(+) T-cell subsets that have divergent transcriptional programs and that display mechanistically distinct cytotoxic potentials. Importantly, an individual's Vδ2 profile predicts defined effector capacities, demonstrated by contrasting mechanisms and efficiencies of killing of a range of tumor cell lines. In short, these data support patient stratification to identify individuals with Vδ2 profiles that have effector mechanisms compatible with tumor killing and suggest that tailored Vδ2-profile-specific activation protocols may maximize the chances of future treatment success.

Separation of functionally divergent muscle precursor cell populations from porcine juvenile muscles by discontinuous Percoll density gradient centrifugation.

BACKGROUND: Satellite cells (SC) and their descendants, muscle precursor cells (MPC), play a key role in postnatal muscle development, regeneration, and plasticity. Several studies have provided evidence that SC and MPC represent a heterogeneous population differing in their biochemical and functional properties. The identification and characterization of functionally divergent SC subpopulations should help to reveal the precise involvement of SC/MPC in these myogenic processes. The aim of the present work was therefore to separate SC subpopulations by using Percoll gradients and to characterize their myogenic marker profiles and their functional properties (adhesion, proliferation, and differentiation). RESULTS: SC/MPC from muscles of 4-day-old piglets were isolated by trypsin digestion and enriched by Percoll density gradient centrifugation. A mixed myogenic cell population was obtained from the 40/70% interface (termed: mixed P40/70) of a 25/40/70% Percoll gradient. Thereafter, by using a more stepped 25/40/50/70% Percoll gradient, mixed P40/70 was divided into subpopulation 40/50 (SP40/50) collected from the 40/50% interface and subpopulation 50/70 (SP50/70) collected from the 50/70% interface. All three isolated populations proliferated and showed a myogenic phenotype characterized by the ability to express myogenic markers (Pax7, MyoD1, Desmin, and MyoG) and to differentiate into myotubes. However, compared with mixed P40/70, SP40/50 and SP50/70 exhibited distinct functional behavior. Growth kinetic curves over 90 h obtained by the xCELLigence system and proliferation assays revealed that SP40/50 and mixed P40/70 constituted a fast adhering and fast proliferating phenotype. In contrast, SP50/70 showed considerably slower adhesion and proliferation. The fast-proliferating SP40/50 showed the highest myogenic differentiation potential with higher fusion rates and the formation of more middle-sized and large myotubes. CONCLUSIONS: The described Percoll density gradient centrifugation represents a useful tool for subdividing pig SC/MPC populations with divergent myogenic functions. The physiological role of SC subpopulations during myogenesis and the interaction of these populations can now be analyzed to a greater extent, shedding light on postnatal growth variations in pigs and probably in other animals.

Novel Homozygous Variants of SLC13A5 Expand the Functional Heterogeneity of a Homogeneous Syndrome of Early Infantile Epileptic Encephalopathy.

BACKGROUND: Early infantile epileptic encephalopathy 25 (EIEE25) is a distinct type of neonatal epileptic encephalopathy caused by autosomal recessive mutations in the SLC13A5 gene. SLC13A5 encodes a transmembrane sodium/citrate cotransporter required for regulating citrate entry into cells. METHODS: Four families with recessively inherited epileptic encephalopathy were sequenced by clinically accredited laboratories using commercially available epilepsy gene panels. Patients were examined by a neurologist and were clinically diagnosed with infantile epileptic encephalopathy. RESULTS: We present four families with global developmental delay, intellectual disability, and defective tooth development with four novel homozygous mutations in SLC13A5. The neurological examination showed spastic quadriplegia with increased deep tendon reflexes. Brain magnetic resonance imaging showed nonspecific signal abnormality of the bilateral hemispheric white matter. Despite similar clinical features, the conditions were based on different molecular mechanisms acting on SLC13A5 (abnormal splicing, large-scale deletions, and tandem-residue insertion). CONCLUSIONS: Our results extend the landscape of autosomal recessive inherited homozygous mutations in SLC13A5 that cause a distinctive syndrome of severe neonatal epileptic encephalopathy. Our observations confirm the homogeneity of epileptic encephalopathy and dental abnormalities as a distinct clinical marker for EIEE25 despite the heterogeneous functional and mutational background.

Modelling motor units in 3D: influence on muscle contraction and joint force via a proof of concept simulation.

Functional heterogeneity is a skeletal muscle's ability to generate diverse force vectors through localised motor unit (MU) recruitment. Existing 3D macroscopic continuum-mechanical finite element (FE) muscle models neglect MU anatomy and recruit muscle volume simultaneously, making them unsuitable for studying functional heterogeneity. Here, we develop a method to incorporate MU anatomy and information in 3D models. Virtual fibres in the muscle are grouped into MUs via a novel "virtual innervation" technique, which can control the units' size, shape, position, and overlap. The discrete MU anatomy is then mapped to the FE mesh via statistical averaging, resulting in a volumetric MU distribution. Mesh dependency is investigated using a 2D idealised model and revealed that the amount of MU overlap is inversely proportional to mesh dependency. Simultaneous recruitment of a MU's volume implies that action potentials (AP) propagate instantaneously. A 3D idealised model is used to verify this assumption, revealing that neglecting AP propagation results in a slightly less-steady force, advanced in time by approximately 20 ms, at the tendons. Lastly, the method is applied to a 3D, anatomically realistic model of the masticatory system to demonstrate the functional heterogeneity of masseter muscles in producing bite force. We found that the MU anatomy significantly affected bite force direction compared to bite force magnitude. MU position was much more efficacious in bringing about bite force changes than MU overlap. These results highlight the relevance of MU anatomy to muscle function and joint force, particularly for muscles with complex neuromuscular architecture.