Molecular modeling of lactoferrin for food and nutraceutical applications: insights from in silico techniques

Abstract Lactoferrin is a protein, primarily found in milk that has attracted the interest of the food industries due to its health properties. Nevertheless, the instability of lactoferrin has limited its commercial application. Recent studies have focused on encapsulation to enhance the stability of lactoferrin. However, the molecular insights underlying the changes of structural properties of lactoferrin and the interaction with protectants remain poorly understood. Computational approaches have proven useful in understanding the structural properties of molecules and the key binding with other constituents. In this review, comprehensive information on the structure and function of lactoferrin and the binding with various molecules for food purposes are reviewed, with a special emphasis on the use of molecular dynamics simulations. The results demonstrate the application of modeling and simulations to determine key residues of lactoferrin responsible for its stability and interactions with other biomolecular components under various conditions, which are also associated with its functional benefits. These have also been extended into the potential creation of enhanced lactoferrin for commercial purposes. This review provides valuable strategies in designing novel nutraceuticals for food science practitioners and those who have interests in acquiring familiarity with the application of computational modeling for food and health purposes.


Introduction
Food loss and food waste, collectively known as food wastage, are global challenges, resulting in the loss of about USD 1 trillion each year (Scialabba 2014).In economically prosperous countries, including Australia, the expenditure involved in confronting this issue is approximately AUD 20 billion annually (Government Australia 2017).Aside from jeopardizing the economy, food wastage also contributes to environmental problems, including 13% of the total greenhouse gas emissions These are primarily as a result of the production of methane gas from livestock, including emissions from cattle and fertilizers (Council Climate 2016;Yuan et al. 2018).In the food industry, the manufacturing of cheese inadvertently yields waste in the form of whey (Lappa et al. 2019).Previous work has focused on salvaging milk solids and refining them into food ingredients, encompassing whey protein concentrate (WPC) and pure lactose powder (Das et al. 2016).In addition to these ingredients, there are also some valuable raw materials, which may be refined into functional foods, which are left in the whey component of milk, including lactoferrin, serum albumin, and immunoglobulin (Patel 2015).
The term functional food was first coined and promoted in the 1980s in Japan in order to mitigate the increasing health-care costs as well as to enhance the health of the aging Japanese population (Henry 2010).According to the International Life Science Institute, functional food is defined as "foods that by virtue of the presence of physiologically active food components provide health benefits beyond basic nutrition" (Crowe and Francis 2013).Lactoferrin (Lf) is a bioactive compound found in secretory fluids, predominantly in milk, saliva, seminal plasma, and tears (Superti 2020).It has various functionalities, including the capacity to encapsulate other small bioactive molecules for a controlled release of these ingredients, deliver or sequester metallic ions, as well as the ability to act as an antimicrobial agent (Kell, Heyden, and Pretorius 2020;Liu et al. 2018).Due to these health benefits, Lf is often formulated in various food related products, including therapeutic drinks, infant formula powders, and fermented beverages (Niaz et al. 2019).In spite of these advantages, the bioavailability of Lf is relatively low due to its susceptibility to digestive enzymes in the gastrointestinal tract.Additionally, exposure to elevated temperatures during food processing and the gastric pH of the human stomach during digestion result in the denaturation of the protein structure, which has limited its commercial application (Sreedhara et al. 2010).Over the years, there has been a significant body of research, which concentrates on antimicrobials; bioactive proteins; encapsulation; functional foods; lactoferrin; whey proteins maximizing the viability of Lf for oral consumption.Typically, this is achieved by coating Lf with protectants, composed of various hydrocolloids, stemming from carbohydrates, proteins, or lipid nanoparticles, which results the formation various complexes that are more stable under differing temperature or pH conditions (Liu et al. 2018).Nevertheless, to date, the fundamental molecular mechanism, mediating its protection by encapsulants under various processing regimes is poorly understood.
Molecular modeling and simulations are powerful techniques, commonly used in the area of materials science, engineering, structural biology, and drug discovery (Aminpour, Montemagno, and Tuszynski 2019;Keil 2018).They have proven to be useful tools in investigating the molecular behavior of certain systems under specific conditions, which is difficult to be detected solely by experimental procedures (van Mil, Boerwinkel, and Waarlo 2013).They are also relatively rapid approaches to predict the binding strength, reaction mechanism, and the structural features of a system, including those in food (Pitsillou et al. 2020;Chen et al. 2019).For instance, in the field of pharmaceutical science, these approaches, verified by an enzyme-linked immunosorbent assay have successfully been implemented to screen and examine the interaction between 300 small molecules and a target enzyme (Pitsillou et al. 2020).These also encompass predictions of binding interactions between antiviral drugs as well as dietary compounds to the main protease involved in the replication of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (Pitsillou et al. 2020).In another example, molecular modeling and simulations could also be potentially used to create enzymes that could be tailored to fit industrial scale production (Chen, Zhang, and Mu 2021).Hence, molecular modeling and simulations can be used as complementary methods to understand the molecular behavior, function, and chemical reactions of numerous systems under various conditions, including food systems.

Scope of this review
In this review, the general structural and functional properties of Lf are elaborated, followed by an overview of molecular modeling and simulations to aid readers with limited or no background in molecular modeling and simulations.Specifically, molecular dynamics (MD) simulations are highlighted as the primary simulation technique used to study Lf and other milk proteins.Direct computational studies involving the encapsulation of Lf and other in silico studies pertaining to stability of Lf for food application are reviewed.These include the changes in structure and function of Lf under food processing regimes.The studies evaluating the interaction of molecules with Lf are also reviewed, with a special emphasis on the binding process.The review concludes with a summary and outlook for future studies.It is hoped that this review will bridge the gap between experimentalist and computational experts, working in the field and will assist in providing future directions for this research area.

Structure and function of lactoferrin
Lf is a multifunctional globular protein, with a molecular weight of ∼80kDa.It is a single chain protein and has two lobes, specifically the N-lobe or the left hemisphere (residue 1 − 333) and the C-lobe or the right hemisphere (residue 345 − 389), which are bound together by a single salt bridge between 334 − 344 (Haridas, Anderson, and Baker 1995;Moore et al. 1997).These salt bridge interactions, composed of 3-turn α-helices act as a control hinge point during the opening and closing of the lobes of Lf, which occur during iron release or binding (Baker and Baker 2009).Each hemisphere can also be compartmentalized further into two sub-units, encompassing N-1 (residue 1 − 90 & 251 − 333) and N-2 (residue 91 − 250) for the left lobe as well as C-1 (residue 345 − 431 & 595 − 676) and C-2 (residue 432 − 592) for the right lobe (Figure 1A).Depending on the mammalian species, there are differences in the amino acid composition of Lf.For example, bovine Lf (bLf) has 689 amino acid residues, whereas human Lf (hLf) has 691amino acid residues (Haridas, Anderson, and Baker 1995;Moore et al. 1997).Despite these subtle differences, sequence analysis showed that the amino acid sequences of bLf and HLf share approximately 68% identity, indicating that the structure and primary biological activity is conserved between bLf and HLf (refer to Figure S1 in Appendix A, Supplementary data) (Yount et al. 2007;Madeira et al. 2019).Nonetheless, these differences will affect the molecular structure of Lf.The superimposed structures, comparing bLf and hLf can be seen in Figure 1B, whereas other genetic sequence alignment comparing Lf, sourced from diverse mammals can be seen in Table 1 (Roberts et al. 2006;Sharma et al. 1999;Karthikeyan et al. 2000;Kumar, Yadav, and Singh 2002;Khan et al. 2001).

Antibacterial properties of lactoferrin
Superimposition of the overall tertiary structures of bLf and hLf indicates that their structures are similar, thereby, supporting the results from the sequence analysis (Roberts et al. 2006).Furthermore, as seen in Figure 1C, embedded within the tertiary structure, there are two peptide segments identified to have antibacterial function.These are known as lactoferricin and lactoferrampin (Bellamy et al. 1992;van Der Kraan et al. 2004).Lactoferricin is the first peptide sequence, isolated in 1991 through pepsin digestion that is known for its potent antibacterial properties against Gram positive and negative bacteria, whereas lactoferrampin is the more recent peptide found in 2004 that has also exhibited both candidacidal and antibacterial activities against various pathogenic bacteria with the exception of those who could ferment for their energy supply (Bellamy et al. 1992;van Der Kraan et al. 2004).
When comparing the antibacterial peptide regions between bLf and hLf, there are some differences in not only the number, but also the composition of the amino acid residues (Figure 1C).In bLf, the segment containing lactoferricin, one of the anti-bacterial peptides in bLf has 25 amino acid residues (residue 17 to 41), whereas in hLf, the number of the amino acid residues is 47 (residue 1 to 47) (Bellamy et al. 1992).Similarly, the number of amino acid residues making up bovine based lactoferrampin, the other anti-bacterial domain in bLf, is 20 (residue 265 to 284), whilst in hLf, there are 17 amino acid residues (residue 269 to 285) (Haney et al. 2009;van Der Kraan et al. 2004).The amino acid profile of both lactoferricin and lactoferrampin derived from humans and bovine can be seen in Table 2.
The molecular mechanism underlying this bactericidal activity involves the positively charged surface of these antibacterial peptides, which are positioned in the N1-lobe of Lf (Figure 1C).Studies showed that this region can hinder the lipopolysaccharide layer of Gram-negative bacteria from collecting magnesium and calcium ions, which are two pivotal minerals for bacteria to grow (González-Chávez, Arévalo-Gallegos, and Rascón-Cruz 2009).A similar trend was also exhibited in the binding of Lf to the cell wall of Gram-positive bacteria, where the positively charged region in Lf binds to the anionic surface of the peptidoglycan layer as well as the teichoic acids (González-Chávez, Arévalo-Gallegos, and Rascón-Cruz 2009).Since bLf has a more positively charged surface than its human counterpart, it is proposed that bLf is more effective than hLf in inhibiting the growth of these pathogenic bacteria (Wang et al. 2019).This is beneficial to the food industry as bovine derived Lf is technically more available in the market than hLf.Overall, the molecular mechanism underlying the antibacterial function of Lf is by disrupting the membrane layer and affecting the homeostasis of these pathogenic bacteria.
In addition to hindering the growth of bacteria, which caused food poisoning, Lf was found to be potent in decreasing the rate of biofilm formation (García-Montoya et al. 2012).For example, in vitro studies have demonstrated the capability of Lf in retarding the production of biofilm by various primary and opportunistic pathogens, including Bacterioides fragilis, B. thetaiotaomicron, Pseudomonas aeruginosa PA01, L. monocytogenes, L. innocua, Streptococcus pneumoniae, and S. aureus (Angulo-Zamudio et al. 2019;de Sá Almeida et al. 2020;Ramamourthy and Vogel 2021;García-Borjas et al. 2021;Shahidi et al. 2020).Despite this functionality, the molecular mechanism responsible for the inhibition of biofilm formation by Lf is not well understood and it is still an on-going field of research (Ammons and Copié 2013).One mechanism proposed for this delayed formation of biofilms is that Lf acts by preventing the adhesion of lectin-dependent bacteria or stimulating the motility of certain species of bacteria (Ammons and Copié 2013).
Practically, this shows a promising application of Lf as a therapeutic agent for tissue engineering in the medical field, including dental, hip, and knee implants where these medical devices are prone to biofilm contamination (Pall and Roman 2020).This application is also relevant for the dairy industry where there have been concerns regarding the formation of biofilms on the surface of stainless-steel pipes, milk silos, and milk tanks (DeFlorio et al. 2021).

Other antimicrobial properties of lactoferrin
As well as being a natural bactericidal agent, Lf can also inhibit the growth of other infectious microorganisms, including viruses and fungi.As an anti-fungal agent, Lf has been proven effective against a broad-spectrum of fungi, including Aspergillus niger, Aspergillus fumigatus, Candida albicans, Candida krusei, Cryptococcus neoformans, and Saccharomyces cerevisiae (Fernandes and Carter 2017;Lahoz et al. 2008;Nikawa et al. 1993;Zarember et al. 2007).This anti-fungal activity is also present in the anti-bacterial peptide regions in Lf.It was found that lactoferricin can curb the growth of various fungi, including C. albicans, C. Krusei, A. versicolor, C. neoformans, S. cerevisiae, Penicillium pinophilum, and P. vermiculatum (Fernandes and Carter 2017;Bellamy et al. 1993Bellamy et al. , 1994)).In contrast, lactoferrampin is only able to delay the growth of C. albicans (Bolscher et al. 2006;Fernandes and Carter 2017).Lactoferricin is thus, more versatile in repressing the growth of various fungi than lactoferrampin.From these results, it could be argued that these peptides are most likely the responsible regions in Lf that perform this inhibitory activity.
It was also found that different forms of Lf have differing efficacy against yeast and molds.Similar to the advantage of apo-Lf in depriving iron from iron dependent bacteria, it was discovered that the iron depleted form of Lf is more potent in decreasing the growth of these iron reliant fungi (Fernandes, Weeks, and Carter 2020).Fernandes, Weeks, and Carter (2020) determined that when Lf is saturated with iron, there is a significant increase in the minimum inhibitory concentration (MIC), which is above the tested concentration (>256 μg/ml).Conversely, the MIC for apo-Lf, sourced from Sigma-Aldrich ranges between 16 to 32 μg/ ml, whereas the MIC spectrum for apo-Lf, obtained from dairy companies is from 128 to 256 μg/ml (Fernandes, Weeks, and Carter 2020).The families of fungi, which apo-Lf could effectively inhibit includes Mucoromycetes, Tremellomycetes, Dothideomycetes, Sordariomycetes, Eurotiomycetes, and Sacchromycetes (Fernandes, Weeks, and Carter 2020).This overall shows that apo-Lf is more functional in countering fungal infection compared to its other forms.
Furthermore, Lf also been demonstrated to be efficacious in preventing the contraction of numerous viral diseases (Wakabayashi et al. 2014).For instance, Lf could obviate the growth of projectile vomiting virus, specifically norovirus in children as well as rotavirus, one of the most common viruses associated with gastroenteritis (Tsukahara et al. 2020;Superti et al. 1997).Lf is also effective against other types of viral infections, including hepatitis c, hepatitis b, alphavirus, hantavirus, human papillomavirus, herpes, poliovirus, and human immunodeficiency virus (van der Strate et al. 2001;Berlutti et al. 2011).This advantage is also extended other respiratory related viruses, including respiratory syncytial virus, influenza, adenovirus, as well as enterovirus (Wakabayashi et al. 2014).The specific regions in which Lf binds to various viruses is contentious and a subject of on-going research.According to a study conducted by Florian et al. (2013), it was found that the peptides derived from the N-lobe of Lf are effective in inhibiting the growth of hepatitis b.They argued that the cationic surface of this region, specifically the GRRRR cluster or the pentapeptide region, exclusive to hLf (residue 1 to 5) (Rybarczyk et al. 2019) interacts with the negatively charged residues of the capsid of this virus and inhibits the growth of hepatitis b by up to 75% (Florian et al. 2013).A similar trend was also observed with the interaction between this domain of Lf with cytomegalovirus as well as hepatitis c virus where there is an electrostatic interaction binding the N terminal of Lf with these viruses (Redwan et al. 2014;Harmsen et al. 1995).Conversely, Scala et al. (2017) found that in the case of the influenza virus, particularly Influenza A virus, the C-lobe of Lf is more effective in limiting the growth of this virus.Specifically, the residues 418 to 429 (SKHSSLDCVLRP) of the C terminus of Lf is able to interact with this virus (Scala et al. 2017).A representation of this segment of Lf is shown in Figure 2A.
Recently, there have been some novel studies on the application of Lf to treat COVID-19, a disease caused by SARS-CoV-2.It was found that the combination of Lf encapsulated with liposomes, vitamin C, and zinc solution could alleviate the severity of the infection as well as reduce the duration of the infection (Serrano et al. 2020).According to Chang, Ng, and Sun (2020), one of the mechanisms on how Lf prevents SARS-CoV-2 from infecting human cells is via direct binding to the spike proteins.Inevitably, this causes difficulty for the virus to bind to the host ACE-2 receptor and seize control of the machinery of the human cells for subsequent replication.The specific molecular mechanism proposed by Miotto et al. (2021) describes that the N lobe and, even more so, the C lobe of Lf as being more available for the binding with the spike protein due to less glycosylation sites being present at these regions of Lf.This result is corroborated by those of Casalino et al. (2020) regarding the significant roles of glycans in modulating the recognition of ACE-2 receptors by SARS-CoV-2.These glycans complementarily sits on the surface of the spike proteins, thereby, creating a high concentration of glycosylation sites (Casalino et al. 2020).As a result of this difference in charge, the C lobe of Lf interacts more effectively than the N lobe of Lf.From the results regarding the interaction of Lf with viruses, there are two key regions, located within the N lobe and C lobe of Lf that could interact efficiently with specific viruses.This is in contrast to the interaction between Lf and those of bacteria and fungi where it is predominantly sustained by the iron chelating capability of Lf as well as the two of the peptide motifs, consisting of lactoferricin and lactoferrampin.

Molecular insights into the iron chelating and delivery functions of lactoferrin
The other functionality that Lf is known for is its capacity to transfer iron.As iron is one of the micronutrients in the world that a substantial proportion of the global population is lacking (World Health Organization 2008), Lf has frequently been viewed as one of the delivery vehicles for this essential mineral (Miller 2013).This ability arises from the high affinity of Lf toward iron as demonstrated by the equilibrium dissociation constant value of approximately 10 −20 M (Baker and Baker 2004).Based on the iron concentration level, there are three distinct types of Lf.In its commercial native state, Lf has about 10-20% iron, whilst in its holo state, the % of iron can range from 76 − 100%, leaving the content in the apo state, consists of approximately 5 to 12% of iron (Wang et al. 2019).At the molecular scale, there are therefore three kinds of conformations, which Lf could adopt.These include: (i) holo-Lf, an iron saturated Lf, with two iron ions bound via coordinate bonding to residue sidechain ligands and a carbonate anion within to each of the lobes in Lf, (ii) apo-Lf, an iron free Lf, without any iron ions embedded in its lobes, and (iii) mono-Lf, a partially iron saturated Lf, with one iron ions linked via a coordinate bond to the residue sidechain ligand and a carbonate anion in either the N or the C lobe of Lf.A graphic representation on the binding of iron to the N and the C lobe of Lf is presented in Figure 2B and C, respectively.
The molecular mechanism by which Lf binds iron directly involves the neighboring amino acid residues in the individual lobes, carbonate ion, and a number of hydrogen bonds (Abdallah and Chahine 2000).In the N lobe, His253, Tyr192, Tyr92, and Asp60 are the primary residues that are covalently linked to the ferric cation (Moore et al. 1997;Wang et al. 2019).Similarly, the ferric cation in the C lobe of Lf is covalently bound to His595, Tyr526, Tyr433, and Asp393 (Moore et al. 1997;Wang et al. 2019).In addition to these residues, there is one carbonate ion in each lobe that interacts with the ferric ion, thereby, stabilizing the ferric cations in these pockets (Abdallah and Chahine 2000).Without the presence of each of these carbonate ions, Lf will lose its binding capacity with the ferric cation (Abdallah and Chahine 2000).This binding process is also supported by the α-helix structure, near the iron binding residues, including residues 121 to 131 in the N lobe and 395 to 407 in the C lobe of Lf (Wang et al. 2019).From these processes, both the N and C lobes of Lf will draw in the ferric cation one by one, starting from the C lobe before moving on to the N lobe (Abdallah and Chahine 2000).Finally, there are also non-covalent interaction that assists in the binding process, including the formation of inter-domain hydrogen bonds, which renders holo or iron saturated Lf in a closed conformation (Abdallah and Chahine 2000).As a result of this closed state, holo-Lf has a more compact structure than its apo-form, resulting in a more robust protein when exposed to thermal or environmental stress (Sreedhara et al. 2010).
The mechanistic action behind the iron release of Lf is different compared to the iron binding of Lf.There are some external factors that could be applied to allow Lf to release iron, including pH and the reduction of the ferric cation to ferrous cation (Baker, Baker, and Kidd 2002;Baker and Baker 2009).The pH of the system could be lowered down to around 4.0 for bLf or below 3.5 for hLf, to effectively protonate the amino acid residues that form covalent bonds with the ligand (Baker and Baker 2004;Rastogi et al. 2016;Wang et al. 2019).This weakens the coordination bond between the metals and lactoferrin; thereby, resulting in the release of the ions (Abdallah and Chahine 2000).This protonation also results in the change in the conformation of Lf, leading to the breakage of the intact hydrogen bond that holds the iron in place as well as allowing other metal chelating substrates in the same medium as Lf to capture the leaking iron ions (Abdallah and Chahine 2000).The other key factor that could induce the apo form of Lf is the reduction of the ferric cation to the ferrous cation via a redox reaction (Baker and Baker 2009; Baker, Baker, and Kidd 2002; Wang et al. 2019).This is fundamentally due to the lower binding affinity of Lf to the ferrous form of iron compared to the ferric form (Baker and Baker 2009;Baker, Baker, and Kidd 2002).Despite these shortcomings, compared to the other proteins in the transferrin family, including serum transferrin and ovotransferrin it was found that Lf is more capable of retaining iron under pH stress (Abdallah and Chahine 2000).In the context of practical usage, this shows Lf as a versatile ingredient compared to the other transferrin proteins.In food applications, Lf could also potentially be added as an additional antimicrobial agent or iron vehicle as most food, excluding high acid food, coming from fruits and fermented products that have low acidity, or even slightly basic pH.

Other application of the iron transferring properties of lactoferrin and general applications of the interactions with other cationic metals
There are numerous advantages of using Lf as an iron transporter.Lf was shown in the previous section to be an iron sequestering agent and inhibits the growth of iron dependent pathogens or spoilers.From a nutritional perspective, Lf has also been shown to effectively deliver iron to those deficient of this mineral, including anemic people and pregnant women (Paesano et al. 2010;Miller 2013).This capacity is also extended to people routinely engaged in physical exertion, such as sportspeople, who tend to lose iron through sweats during exercise (Özer and Kirmaci 2010;Tang et al. 2016).In addition, Lf has been demonstrated to have healthy effects by modulating the amount of iron in the human body (Carmona et al. 2017).Although iron is an essential nutrient for cell growth and development, including energy metabolism, excessive iron intake can be toxic for human body, resulting in the production of oxygen radicals (Rosa et al. 2017).Particularly, the free form of iron, also known as the ferrous cation is more hazardous than its ferric cation counterpart because it catalyzes the creation of oxygen radicals via the Haber-Weiss and Fenton reactions (Lepanto et al. 2019;Bokare and Choi 2014).The accumulation of these radicals will not only cause cell apoptosis through damaging cellular components, including the proteins, DNA, and membrane layer, but also ferroptosis where the death of these cells is triggered by an auto-immune response mediated by the glutathione peroxidase 4 (Sun et al. 2020;Lepanto et al. 2019).Eventually, these will result in inflammation.
In practical settings, this iron modulation activity of Lf has been reported to have beneficial effects in retarding the symptoms of neurological illnesses, including Alzheimer's and Parkinson's diseases (Liu et al. 2020;Akilo et al. 2018).Studies have hypothesized and suggested that the pathogenesis of these sicknesses is due to the dysregulation and excess iron uptake in the human brain (Akilo et al. 2018;Liu et al. 2020).Accordingly, Lf could potentially sequester surplus iron, resulting in the reduction of the oxidative stress imparted on the human brain.
Recently, the chelating function of Lf has also been investigated to study the effects of SARS-CoV-2 on human cells (Habib et al. 2021).It was found that the pathogenesis of COVID-19 is also related to elevated iron levels, which indirectly causes inflammation in the human body (Habib et al. 2021).This mechanism seems to be caused by viral entry into human cells, resulting in potential and extended apoptosis of the red blood cells, leading to the dissociation of the porphyrins from hemoglobin and resultant release of iron (Habib et al. 2021;Liu and Li 2020).This leads to increase in free iron content and result in a chain reaction where the amount of ferritin is increased to compensate for the loss of hemoglobin in binding oxygen (Habib et al. 2021;Goldberg et al. 2020;Zhou et al. 2020).Accordingly, the presence of Lf could act as an adjuvant to control this exponential rise in iron through direct binding to iron, and/ or indirect modulation of the cytokines (Habib et al. 2021).In general, whilst undoubtedly holo-Lf is a more robust and suitable ingredient to be used for nutritional purposes compared to its iron depleted form, apo-Lf, the latter form has however been demonstrated to have much more diverse and enhanced functions as a natural anti-microbial agent and therapeutic substance than iron saturated Lf.
Finally, the interaction of Lf could also be extended to other cationic minerals (Honarparvar, Kanchi, and Bisetty 2019).One simulation study demonstrated that Ti 4+ has the strongest binding affinity to Lf (Honarparvar, Kanchi, and Bisetty 2019).This was followed by Fe 3+ , Cr 3+ , Co 2+ , and Cu 2+ (Honarparvar, Kanchi, and Bisetty 2019).Nevertheless, recent experimental studies also highlighted that Lf can also bind to Zn 2+ , as shown by the exothermic reaction (Tang and Skibsted 2016).Practically, a composite biomaterial, composed of Lf and titanium could be used to as a next generation dental implant with antimicrobial properties (Yin et al. 2020).Furthermore, as Lf has a relatively binding affinity against Cr 3+ , it was discovered that it could potentially be used as an immunomodulatory therapy for chromium-induced kidney injury (Hegazy et al. 2016).Interestingly, it was also found by supplementing Lf with Cu 2+ and Mn 2+ , the anticancer properties of Lf is elevated, particularly in adenocarcinoma gastric cell lines (Bo, Li, and Zhao 2019).Overall, the findings indicated that depending on the attached cationic minerals, the properties of Lf could be tailored to fit various purposes, including health and antimicrobial capacities.

Overview of molecular modeling and simulations
There are many simulation techniques, which can be classified based upon the size of the system, time scale of simulation runs, accuracy, computational resources required, as well as the complexity of the calculation (Figure 3).
From the diagram, finite elements methods, also known as computational fluid dynamics can be applied as one of the simulation techniques used to model large systems exceeding 10 µm.For example, these calculations could be used to model heat flow in ovens and chillers, flow of food in continuous-flow systems (Scott and Richardson 1997).In contrast, Monte Carlo based simulations are used to model systems with a size between 100 nm − 10 µm.This type of simulation is particularly useful to monitor the growth as well as calculate the thermodynamic properties of structures with various crystalline forms (Brito and Cândido 2020).It has also been demonstrated that in conjunction with atomistic models, it could be used to model biological systems, particularly proteins and their binding with small molecules, which may be useful to model certain proteins in food (Cabeza de Vaca et al. 2018).Nevertheless, in general, these two methods are primarily useful in the field of materials science and engineering.
In order to obtain an in-depth molecular understanding of biological systems, specifically in food research related areas, there are three common types of simulation procedures, which could be considered.These are coarse-grained (CG) molecular dynamics simulations, atomistic molecular dynamics (MD) simulations, and quantum mechanics (QM) based simulations (highlighted in the black box of Figure 3).The CG technique, such as the approach adopted by the MARTINI forcefield could be used to model various food macromolecules, including carbohydrates, proteins, and lipids (Kmiecik et al. 2016).Recent work utilizing the MARTINI forcefield have shown promising application in modeling the lipid polysaccharide layer of Gram-negative bacteria and the phospholipid bilayer in the cell membrane, which are particularly useful in the field of food microbiology (Li, Fang, et al. 2021;Shearer et al. 2020).Moreover, in the context of proteins, CG models could also be used to understand the association among different types of proteins, including the dimerization or self-assembly of proteins in membranes, which are useful in designing novel therapeutical agents (Lamprakis et al. 2021).In contrast, QM calculations are used to model the interaction of miniscule systems, especially those involving covalent interactions (Levine and Head-Gordon 2020).For food applications, this type of calculation could potentially be used to elucidate the interaction of various food molecules, including many food proteins, such as those from milk, eggs, and soy, which can be cross-linked with polyphenols and phenolic acids (Keppler, Schwarz, and van der Goot 2020).Hence, QM calculations could give insights regarding the molecular mechanism behind the binding with small molecules, which results in increased functionality of these proteins, including enhanced emulsifying capacity as well as controlled-release of the dietary compounds (Keppler, Schwarz, and van der Goot 2020).
Despite the advantages of CG models and QM calculations, there are some disadvantages, which may limit the application of these simulation techniques for food purposes.In the food industry, thermal or non-thermal processing is commonly implemented to ensure the safety of food (Zhang et al. 2020).Inevitably this leads to changes in the molecular behavior of the food components, particularly proteins as they are prone to undergo denaturation under non-physiological conditions (Narayan, Bhattacharjee, and Naganathan 2019).Unfortunately, many CG models could not effectively model the effects of environmental changes, including thermal processing, on the intramolecular and structural changes because they simplify the representation of the macromolecular system, which replace all atomistic details of the studied systems with larger interaction sites and often require artificial force restraints to ensure molecular structural integrity (Singh and Li 2019).Thus, quantitative thermodynamic properties and the molecular details of the binding process of the system being modeled cannot be fully determined using most CG approaches (Oprzeska-Zingrebe and Smiatek 2019).Although QM calculations could be used to accurately model systems, involving covalent interactions, the computational expense involved and the size of the system often restrict the capability to employ them for large systems, including proteins (Morawietz and Artrith 2021).One well-known proprietary software used to model systems based upon the principle of QM calculations, called Vienna ab initio simulation package (VASP) could only model up to thousands of atoms, rendering it unfeasible to model the complete structure of food macromolecules, especially protein-based matrices (Kresse andFurthmüller 1996a, 1996b;Kresse and Hafner 1993).Conversely, with increasingly more powerful high performance computing (HPC) technology, MD simulations offer promising application in food as it can model larger molecular systems, incorporating more macro-components, including carbohydrates, proteins, enzymes, and lipids (Li et al. 2020;Shukla and Tripathi 2020).For instance, recent studies have demonstrated the application of MD simulations to understand the molecular behavior of globulin and prolamin, which are two storage proteins, commonly found in plants, including cereal grains (Fang et al. 2022;Hajjari and Sharif 2021).In essence, this technique compared to other simulation techniques could provide valuable insights into understanding the structural and functional changes of the food components, including the physiochemical properties for engineering novel food materials with tailored functionalities (Chen et al. 2019;Chen, Zhang, and Mu 2021).Hence, in this review, the atomistic MD simulation method is highlighted as the primary method to understand the structural changes in food molecules.This is a potentially cost-effective and novel technique in the field of food science and technology and its application to evaluate the stability of Lf under diverse food processing conditions will be reviewed.In addition, the capacity of atomistic MD simulations to elucidate molecular-level details in food molecules enable such techniques to complement higher-level simulations and models commonly employed in food science.For example, atomic simulations can providing the microscopic parameters required for larger-scale continuum or mathematical models which simulate macroscopic physical processes occurring in food systems, such as mass and heat transfer, texture, and rheological properties (for a recent review of larger scale modeling in food science and engineering, see Vitrac, Nguyen, and Hayert (2022)).These macroscopic properties are vitally underpinned by microscopic, molecular-level processes which can only be revealed by atomistic MD simulations.

Theory of molecular dynamics simulations
MD simulation, also known as the computational microscope, is a technique stemming from the principles of classical mechanics (Sekercioglu and Duysak 2009).In a molecule, each atom is treated as a sphere, whereas the bonds connecting different atoms are regarded as springs (Lamberti et al. 2002).A schematic diagram showing an overview of MD simulations is presented in Figure 4.
There are three different sets of data required to run MD simulations for proteins (Sonavane et al. 2014).Firstly, a file, typically in the pdb format, is essential to depict the initial 3 D coordinate of the molecule.These starting structures are generally obtained from experimental data, specifically XRD, Cryo-EM, or NMR, which are accessible from the Protein Data Bank (PDB) (https://www.rcsb.org/).In addition, many biomolecular MD simulation packages make use of a topology file, which define the bonding geometry of the molecule and the force field parameters to be used.These describe the inter-atomic potentials, partial charges, and other nonstandard geometrical restraints, including distance restraints between non-bonded atoms.In general, parameters for a number of commonly-used force fields are obtained from ab initio calculations and/or experimental data (Fröhlking et al. 2020).Many commercially available software programs have been used to model the topology of biomolecules, including CHARMM, GROMOS, and AMBER (Guvench and MacKerell 2008).Each program typically uses specific sets of force field, with its own parameterization.However, these programs can generally produce similar trends of the modeled systems (Hospital et al. 2015;Rueda et al. 2007).The equation for a molecular forcefield can be seen in Equation 1: where E bonded refers to the total energy potential that contributes from bonded interaction, E nonbonded is from non-bonded interaction, leaving E other , which are additional potential energy functions used in specific force fields to augment the accuracy of the calculation (Mackerell Jr. 2004).
The E bonded equation is expanded regarding the specific bonded interaction potentials and can be seen in Equation 2: The first term in Equation 2 is analogous to the potential energy of springs, originating from Hooke's law, which describes the stretching of the bonds using a quadratic equation (refer to FigureS2 in Appendix A, Supplementary data).The K i are spring constant terms, with i and i 0 , which represent the actual bond length and the equilibrium bond length.The second term, also stemming from Hooke's law constrain the bending angle (θ ) of two bonds for three consecutively bonded atoms, with K θ and θ 0 serving as the spring constant and equilibrium bond angle, respectively (refer to FigureS3 in Appendix A, Supplementary data).The third term in the mathematical expression takes account of the dihedral angle, formed over four covalently linked atoms (refer to FigureS4 in Appendix A, Supplementary data).A trigonometric cosine formula is used to model the periodicity of the dihedral rotation, with θ 0 functioning as the constant and n serving as the periodicity.The ∅ is the actual angle, whereas the ∅ 0 is the reference angle.Finally, the last term in E bonded is used to maintain the configuration of the chiral as well as to designate the planarity of specific bonded atoms, including functional groups involving carbon, oxygen nitrogen atoms with a hybridization of sp 2 , as well as those of carbonyl groups and aromatic rings (refer to FigureS5 in Appendix A, Supplementary data).This mathematical term ensures the accuracy of the calculation because the dihedral angle equation is not able to sufficiently take account the planarity of sp 2 geometry (Leach 2001).Similarly, K ω is the spring constant, while ω is the corresponding deviation, with ω 0 as the reference point.
The equation describing non-bonded interactions, which is vital to modeling the inter-molecular interactions, is illustrated in Equation 3: (3) Two of the most common non-bonded interactions include van der Waals and electrostatic interactions (Guvench and MacKerell 2008).The first mathematical expression refers to the van der Waals interaction, approximated using the Lennard-Jones potential.ε ij is the coefficient of the Lennard-Jones potential, with σ ij as the distance where the Lennard-Jones potential is at the minimum point and r ij as the inter-atomic distance (refer to FigureS6 in Appendix A, Supplementary data).The second mathematical term includes the electrostatic interaction based upon Coulomb's law (refer to FigureS7 in Appendix A, Supplementary data).The charges between particles are defined as q i and q j , respectively, with ε 0 being the Coulomb's electric constant.
The E other term can include the Urey-Bradley angle term, exclusive in the CHARMM force field which introduces additional degrees of freedom in order to generate accurate vibration spectra of the molecule (Vanommeslaeghe et al. 2010).This equation is depicted as Equation 4.
In Equation 4, the K UB is the spring constant, whereas the s and s 0 are described as the distance and the reference distance between the atoms, respectively (refer to FigureS8 in Appendix A, Supplementary data).
Other terms that are usually added into force fields in order to fit to the molecular vibration spectra during parametrisation include hydrogen bonding potential and the polarization energy (Gavezzotti 2006;Warshel, Kato, and Pisliakov 2007).These simplified equations are shown in Equations 5 and 6, separately.
In a hydrogen bonding potential, A and B, refer to the two particles, usually between F, O, N, and the hydrogen itself (refer to FigureS9 in Appendix A, Supplementary data).The r AD is the distance between the donor and acceptor, whereas the θ AHD is the angle between the donor and the acceptor.
The polarization energy is designed to enable classical force fields to model charge redistribution in atoms and molecules when exposed to different electrostatic environments, which induces temporary dipoles responsible for London dispersion or van der Waals interactions (Ahlström et al. 1989;Rick and Berne 1996;Warshel, Kato, and Pisliakov 2007).In Equation 6, E PE describes the magnitude of the induced point dipole of a polarizable atom, whilst α i and E i are defined as the polarization of the atom and the electric field of an atom i, respectively.Following from the above mathematical terms regarding the molecular mechanics model employed in classical MD simulations, the physical movement of these atoms can be calculated by applying the Newton's second law of motion, as shown in Equation 7.
By integrating Equation 7 as well as calculating the acceleration of each atom, the trajectory of the whole system is updated periodically in order to adjust the position of each atom over the discretized simulation time (Leach 2001).The equation of this integrated formula is shown in Equation 8.
where F i is the force, acting on an atom i , r i is the position, and m i is the mass of atom i.To expand this equation further, given that not only the total force for all atoms at a time t is computed as the vector sum of its interactions with the other atoms, but also their specific position and the velocity at a particular time t , the updated position and velocity are then recalculated at t t G , in which δt is the interval of time between two timesteps of the MD simulation.During this calculation, the force is also assumed to be constant, leading to new positions and velocities at time t t 2G (Leach 2001).In addition, to minimize the instability of the system, this time step must be a magnitude smaller than the fastest time scale motions within the system, which are typically vibrational motions of bonded atom.However, it must not be too small to prevent any unnecessary computational costs (Hospital et al. 2015;Kim 2014).A time step between 1 and 2 fs is generally selected for running MD simulations (Leach 2001).
In order to determine new positions and velocities at t t G , integration algorithms are implemented.Among these vast, widely available protocols, the Verlet algorithm and the leap-frog algorithm are two of the most commonly used methods to integrate Newton's equation (Leimkuhler, Reich, and Skeel 1996;Chen 2018).The Verlet algorithm (Equation 9) computes the new position of each atom by subtracting the value of the past position [ r t t G ] from that of two times the value of the present position [ r t ].A polynomial expansion, in the form of Taylor's series is often used to complement the Verlet algorithm to calculate the velocity of the system.This calculation can be achieved by subtracting the r t t G from the r t t G before being divided by 2δt .
The benefit of this approach lies on its simplicity as well as improving the computational demand.However, this method has limited robustness and an implicit velocity term, rendering it difficult to obtain the value of the velocity without prior known information of the before and after positions (Leimkuhler, Reich, and Skeel 1996;Chen 2018).To overcome this shortcoming, the leap-frog algorithm was developed, based upon the following relationship between Equation 11 and Equation 12: As seen in the above equations, the velocity with the acceleration at present time point >Gt t a ] .
Subsequently, the position r t t G is obtained from the calculated velocity, with r t as the initial position.To acquire the full-step velocity, Equation 13 may be used: For this equation, the velocity is calculated at t 1 2 Gt , although the result obtained for the position is at [G ta ( )] t , thereby, earning the name leap-frog algorithm.One primary advantage of this approach is that all the parameters concerning the velocity are directly calculated.Nonetheless, the position could not be obtained from the same time as the velocity of the system (Chen 2018).Overall, each algorithm has its own advantages and disadvantages (Leach 2001).Depending on the computational cost and the precision of the data required, some algorithms are preferred over others (Leach 2001).Other algorithms, currently being developed and used in certain simulations include the predictor-corrector protocols, velocity Verlet approach, and the Runge-Kutta methods (Leach 2001;Leimkuhler, Reich, and Skeel 1996;Chen 2018).
The aforementioned mathematical expressions demonstrate the fundamental theories and algorithms behind MD simulations.The heuristic approach on how to perform MD simulations can be seen in Figure 5.

Workflow of molecular dynamics simulations
Practically, the application of MD simulations in biomolecular science often involves the use of MD simulation software packages.There are many software packages that could be used to model biological macromolecules, including Nanoscale Molecular Dynamics (NAMD) and Assisted Model Building with Energy Refinement (AMBER) (Phillips et al. 2020;Case et al. 2005).Nonetheless, one commonly used software is the GROningen Machine for Chemical Simulations (GROMACS), a state-of-the-art, open-source, and community supported set of codes capable of modeling protein, lipid, and carbohydrate molecules (Abraham et al. 2015;Hess et al. 2008).Prior to running GROMACS, the molecular structure of interest ought to be checked to ensure that there are not any missing residues.This can be achieved by using Visual Molecular Dynamics (VMD), Pymol, or many other visualizing software designed specifically to evaluate the structural properties of diverse molecules (DeLano 2002;Humphrey, Dalke, and Schulten 1996).Having confirmed that the structure initially acquired from experiments is intact, the topology of the molecule could be generated via GROMACS, followed by the establishment of the periodic boundary conditions (PBC) to a represent bulk system using a relatively small number of particles (Leach 2001).The size of the periodic box must also be adequate to prevent any mirror-image structures from interacting with one another in the x,y, and z directions that could potentially lead to artifacts (Reed and Flurchick 1996).The term PBC can be interpreted as placing the molecule in a self-contained 3 D supercell, surrounded by 26 mirror-images (refer to FigureS10 in Appendix A, Supplementary data).Different shapes of periodic supercells may be employed, including cuboid, rhombic dodecahedron, and truncated octahedron.Depending on the molecules and the computational efficiency desired, some periodic box geometries are preferred over others (Wassenaar and Mark 2006).In general, the cubic cell is the simplest type of PBC to be used (Leach 2001).
After placing the molecule in a supercell, water and salt are usually added to fill the vacuum space inside the cell, neutralize the system, and provide a salt concentration similar to the physiological or experimental conditions modeled.There are various types of water models that could used to solvate the system, based on the force fields.For example, the TIP3P water model is tailored for both CHARMM and AMBER force fields, whereas the SPC is suitable for the GROMOS force field (Berendsen et al. 1981;Jorgensen et al. 1983).In addition to the compatibility of each water model with the specific force field, there are some differences between the two water models.For instance, the geometry hydrogen charges within the Lennard-Jones parameters of SPC and TIP3P water models are distinct to one another.While the SPC model assumes an ideal values of the O-H bond length of 1 Å and tetrahedral angle of 109.47°centered at the oxygen, the TIP3P model more accurately employs the observed bond length of 0.9572 Å and tetrahedral angle of 104.5° (Leach 2001).These subtle nonetheless important differences will affect the solute and solvent interactions Before running the MD simulations, some preliminary simulations are often implemented to stabilize the system.Energy minimization using the steepest-descent or conjugate gradients algorithms is often used to optimize the geometry of the molecules.This is followed by equilibration simulations, usually performed under constant particle number, volume, and temperature (NVT) conditions using the Berendsen thermostat, followed by simulations under constant particle number, pressure and temperature (NPT) conditions using the Berendsen or Parrinello-Rahman barostats (Berendsen et al. 1984;Parrinello and Rahman 1980, 1981, 1982).These are applied to ensure that the whole system is free from artificial forces due to the semi-crystalline nature of the solvent in the initial system, but also fit to the experimental condition, including allowing for the solvent to form a fluid phase, and for proper modeling of the solvent or environmental pressure as well as controlled temperatures.Finally, having completed these procedures, 'production' MD simulations could be performed, followed by simulation trajectory analyses, employing GROMACS analysis tools as well as VMD, Pymol, or or any of a vast and growing array of structural analysis codes (Abraham et al. 2015;DeLano 2002;Hess et al. 2008;Humphrey, Dalke, and Schulten 1996).

General considerations of other computational methods for food applications
There are other complementary in silico techniques that can be applied in conjunction with MD simulations for modeling Lf in the context of food (Table 3).Commonly used  1. the addition of α-lactalbumin and β-lactoglobulin results in structural cohesion and rigidity.2. the secondary structure of apo-lactoferrin is sustained under simulated freeze-drying conditions.3. the anti-bacterial motifs in apo-lactoferrin is protected under simulated freeze-drying conditions.4. the changes of the distribution of amino acid residues may result in the ease of dispersion in water following the rehydration process. 5. e lectrostatic interaction is unaffected by the freeze-drying temperatures.
( Molecular docking and Md simulations 1. Gln87, Gln89, tyr192, val250, lys301, and ser303 form hydrogen bonding with lactoferrin 2. tyr82, tyr92, and tyr192 interact hydrophobically with lactoferrin 3. Gly83, thr84, Pro88, thr90, His91, arg121, lys210, Pro251, ser252, His253, and arg280 form a van der waal interaction with lactoferrin 4. ser212 forms unfavorable interactions with lactoferrin 5. structural analyses, consisting of root mean square deviation, root mean square fluctuation, radius of gyration, and solvent accessible surface area indicate that the lactoferrin forms a stable complex with theaflavin monogallate (Khan et al. 2021) Free-standing bovine lactoferrin with procyanidin room temperature (25 °C); neutral pH with 10 mM ionic strength Molecular docking 1. Procyanidin was bound to the inner n lobe of lactoferrin.2. Hydrogen bonding and hydrophobic interactions are the primary mode of binding 3. Phe190, Glu80, arg121, Cys9, and Glu15 form hydrogen bond with lactoferrin.4. Phe190,ile11,ala42,val57,lys301,and  membrane protein and the chimeric peptide, the binding loop of Cd8 was not the binding site 3. Hydrogen bonding is the predominant interaction, as verified by previous in vitro study 4. the chimeric peptide binds to the β-sandwich and front layer of the e2 protein 5. β-sandwich domain is the interaction point between the Cd8 membrane protein and e2 glycoprotein 6. Primarily, positively charged and polar amino acid residues of the chimeric peptide interact with the e2 glycoprotein approaches include molecular docking, homology modeling, and de novo protein modeling.Molecular docking software, including Autodock, HEX, and Schrödinger suite are analytical tools that could be used to predict the binding site between protein and ligand, protein and peptide, as well as protein and protein (Ritchie and Kemp 2000;Schrödinger 2020;Trott and Olson 2010).Homology modeling and de novo protein modeling are two methods that are used to predict the structures of proteins, which are not experimentally available (Kuhlman and Bradley 2019).The difference between these two procedures lies in their ways of predicting the structure of proteins.Homology modeling uses a template from a closely related structure, whereas de novo protein modeling is usually used to construct the protein structure without prior knowledge of any closely related proteins (Bordoli et al. 2009;Greener, Kandathil, and Jones 2019).A relevant example from dairy science is the casein proteins in milk.The reliability and stereochemistry of these predicted models are usually analyzed by model evaluation protocols, including PROCHECK and WHATCHECK (Laskowski et al. 1993;Hooft et al. 1996).Some of the freely available servers comprise AlphaFold, SWISS-MODEL, Modeler suite, I-TASSER, and PEP-fold (Waterhouse et al. 2018;Webb and Sali 2016;Maupetit, Derreumaux, and Tufféry 2010;Harpreet, Aarti, and Raghava 2007;Roy, Kucukural, and Zhang 2010;Jumper et al. 2021).Overall, these are the generic methods that could be used to complement MD simulations by producing initial structures that serve as inputs.
To understand how proteins bind to their targeted molecules, free binding calculations could be implemented.These include molecular mechanics Poisson-Boltzmann surface area (MM-PBSA), steered MD simulations coupled with umbrella sampling calculations, linear energy interaction calculation (Baker et al. 2001;Kumari, Kumar, and Lynn 2014;Abraham et al. 2015;Hess et al. 2008) and many other approximate approaches.Depending on the system being studied and the computational cost, some methods are preferred.For instance, MM-PBSA could be used to elucidate the contribution of each residue as well as estimate the energy contribution from electrostatic, van der Waals, polar solvation, and non-polar solvation interactions.Although linear energy interaction calculations are slightly more accurate than the former as it includes induced-fit effects and mobility of water molecules, it could not decompose the energy binding into individual contributions from specific residues of the protein (Rifai et al. 2019;Gutiérrez-de-Terán and Åqvist 2012).Steered MD simulation coupled with umbrella sampling calculation, in contrast, is considered the most accurate method to determine the Gibbs free energy calculation because it involves simulations of the intermediate stages during the binding process, which naturally enables enhancements in the sampling of the structural changes, particularly in high energy regions of the binding pathway that are typically challenging to be accessed using equilibrium MD simulations (Adcock and McCammon 2006;Fogolari, Corazza, and Esposito 2018;Kästner 2011).Nevertheless, this method is more computationally expensive than the other two protocols.Other ways to determine how the complex forms include density functional theory (DFT) calculations, which are based upon the principles of QM to understand the covalent interaction (Kedziora et al. 2016).From a structural perspective, inter-residue contacts between macromolecules can be visualized using network analysis programs, including Gephi and Cytoscape, which translate complex interaction data into an intuitive form and give qualitative insights on the association of residues within the protein matrix and the interacting residues with small molecules (Bastian, Heymann, and Jacomy 2009;Shannon et al. 2003).
Recently, there have been studies investigating the structural dynamics and the encapsulation mechanism of Lf, along with its derivative forms using MD simulations and the above supplementary computational techniques for food and nutraceutical applications.The results of these studies have been summarized in Table 3 and are briefly articulated in the following section.

Application of molecular modeling to probe the stability and encapsulation of lactoferrin under processing conditions
In general, computational approaches have increasingly been established to study the conformation of Lf (refer to Table 3 for further details).For example, the modulation of iron binding as well as the role of carbonate and positively charged residues in controlling this mechanism have been thoroughly elucidated by Anghel, Radulescu, and Erhan (2018).Furthermore, the effects of different solvents, including chloroform and methanol, along with the importance of salt on the structural stability of the anti-bacterial motifs of Lf, specifically lactoferricin have also been investigated by Daidone et al. (2011); Fornili, Pizzi, and Rebeccani (2010).There have also been studies performed in the context of food science.These studies explore the differing effects of food processing regimes, including pasteurization, ultra-high temperature (UHT), and spray drying (Nhan, Rix, et al. 2019).Other conditions, including extreme pH and those encountered in the gastric pH on the stability and the potential functionality of apo-Lf have also been evaluated using MD simulations (Nhan, Small, et al. 2019).Both studies have elucidated fine-structure details, otherwise difficult to probe experimentally, including the differences in conformational response to extreme temperature and pH effects between the N-and C-lobes of apo-Lf (Nhan, Rix, et al. 2019;Nhan, Small, et al. 2019).Moreover, through the use of various computational techniques, these studies have been extended by means of encapsulation of apo-Lf with β-lactoglobulin and α-lactalbumin in order to enhance the stability of apo-Lf when exposed to spray drying (Darmawan et al. 2020), freeze-drying (Darmawan et al. 2021b), and gastric pH (Darmawan et al. 2021c).Studies which investigated the use of Lf as a novel encapsulating agent to have also been investigated.Specifically, these include the potential encapsulation of procyanidin, theaflavin monogallate, naringenin, and LAB with Lf (Li, Dai, et al. 2021;Nunes et al. 2020;Darmawan et al. 2021a;Khan et al. 2021).
As a possible component of novel functional food, there have been studies regarding the use of Lf that promotes its bio-functionality.For instance, Pirkhezranian et al. (2020) have investigated the interaction between the antimicrobial peptides of Lf and bacterial B-DNA.Similarly, Tahmoorespur et al. (2020); (Tanhaeian et al. 2020) have also elucidated the interaction and potential binding capability of a chimeric peptide of Lf with a variety of targeted molecules: the E2 glycoprotein of hepatitis c virus, CD8 human membrane protein, a lipid polysaccharide layer (LPS), and a B-DNA model in bacteria.The N-acylated form of an amidated hexapeptide, RRWQWR-NH2, derived from bovine lactoferricin, and its interaction with the lipid bilayer of bacteria have also been investigated (Romo et al. 2011).Similar membrane binding studies have been extended to the other anti-microbial peptides of Lf, including lactoferrampin (Tsutsumi et al. 2012).In regard to the differing isoforms of lactoferricin, Zhou, Tieleman, and Vogel (2004) have thoroughly investigated and compared the stability of this peptide in its α and β forms.There are also a number of studies, which evaluated the combination between Lf and silver as ions and nanoparticles, highlighting the synergistic effects between these two components as a novel nutraceutical (Nayak et al. 2019;Pomastowski et al. 2016).Overall, it is evident that there is an increasing trend in the use of computational methods, which can complement conventional experimental procedures to understand the molecular stability and binding mechanism of Lf with various molecules.

Conclusions and future directions
The interest in food research that can not only promote technological, but also biological function is clearly evident by the number of studies on functional foods and their health benefits.Whilst this review only focusses on studies involving Lf as a bioactive compound in milk, it is clear that a deep understanding of this protein's numerous biological functions, stability under various conditions, and with other biomolecules is emerging at the molecular and residue level, which can provide insights in the development of novel nutraceuticals.
Overall, computational techniques have proven their advantages in elucidating the interaction between molecules under various conditions.These relatively novel approaches in the field of food science and technology can complement the traditional experimental methods as well as evaluate the stability of compounds when processed in food processing regimes.In this review, the stability and functionality of Lf as well as its interaction with molecules are described in a detailed manner to establish current understanding at the molecular level.Such findings could be beneficial to experimental studies on Lf in explaining the key residues and regions in Lf that have been demonstrated to have functional benefits.For example, modeling has elaborated the role of carbonate ion and Arg121 in modulating the iron delivery capability in Lf.Other modeling studies provided an in-depth molecular level of understanding regarding the mechanisms of the action of specific anti-microbial peptides of Lf, including lactoferricin and lactoferrampin.Studies investigating the effect of pasteurization and ultra-high temperature processing on the structural stability as well as strategies, involving encapsulation to enhance the health benefits of Lf by means of combining with LAB, procyanidin, and naringenin, have also been conducted.Some of these simulated encapsulation studies are also tailored to fit industrial processing considerations under spray drying and freeze-drying as well as gastric digestion conditions encountered as part of human consumption.In the context of evaluating the role of Lf as a novel nutraceutical, recent studies have evaluated the interaction of Lf, including its derivative peptides with a DNA of bacteria, glycoproteins in viruses, human membrane proteins, as well as the lipid bilayers.The work produced from these studies may benefit the pharmaceutical industries, as it highlights the potential use of Lf as natural drugs.In a broader application context, the interaction between Lf and in the form of nanoparticles and cations have also been elucidated.These types of studies may be of benefit for the development of novel materials with antimicrobial properties.Since antibiotic resistance is an emerging health crisis due to increasing robust strains of pathogenic bacteria, alternative treatments, including the use of Lf could potentially mitigate this issue.In general, there is a wide range of the applications of Lf that can benefit food, nutraceutical, and pharmaceutical industries.
Whilst there is a substantial body of computational studies performed on investigating the stability of Lf, there are other various antimicrobial proteins that could be evaluated, including lactoperoxidase and ovotransferrin.Extensive in silico work that could also be performed include studies that examine the interaction of Lf with other ions, including Ti 4+ , Cr 3+ , and Zn 2+ .Nevertheless, it should be noted that there are some limitations that should be considered when running MD simulations, including insufficient samplings of the thermodynamics and kinetic properties of biomolecules.This may limit the capacity of MD simulations to be effectively used to observe dissociation or self-assembly of proteins as well as protein folding as it usually occurs at millisecond timescales, whilst the standard MD simulations are usually performed nanosecond timescales.To compensate for this drawback, often, heavy computational times are required to observe these structural changes.Other ways of minimizing this disadvantage is via applying other advanced MD simulations, including replica exchange molecular dynamics or umbrella sampling calculations, where bias forces are applied in order to accelerate macromolecular changes to allow them to occur on a shorter timescale.
a relatively rapid and cost-effective screening procedure, employing computational approaches may be performed to evaluate whether these compounds are viable to food processing regimes and whether additional protective strategies are required to enhance them, or at least preserve their function.In summary, this review has demonstrated the versatility of computational methods as relatively novel tools to be applied for food and health.Whilst a range of computational techniques have been applied to study Lf, there are other types of functional ingredients that could be further explored using these approaches.

Figure 1 .
Figure 1.(a) structure of lactoferrin with two metal iron.(B) superimposed structures of human (pink) and bovine (blue) lactoferrin.(C) surface charged of lactoferrin at pH 7 and the corresponding lactoferricin (orange & blue) and lactoferrampin (purple & pink).

Figure 2 .
Figure 2. (a) anti-influenza motif in lactoferrin.Cartoon representations of the Fe 3+ (gray) and the Co 3 2-(blue & pink) binding site in the n (B) and C lobe (C) of lactoferrin.

Figure 3 .
Figure 3. various kinds of simulation techniques.

Figure 5 .
Figure 5. Practical workflow of running molecular dynamics simulations.

Table 1 .
structural comparison between lactoferrin sourced from different mammals and lactoferrin obtained from human and bovine.

Table 2 .
Comparison of the amino acid composition for lactoferricin and lactoferrampin sourced from bovine and human lactoferrin.

Table 3 .
independent computational or with experimental studies on lactoferrin for food and nutraceutical application.