Cholesterol affected dynamics of lipids in tailor-made vesicles by ArcVes software during multi micro second coarse grained molecular dynamics simulations
In biotechnology, as well as in nanomedicine, vesicular structures are essential features of the cellular life cycle. Lipid compositions define the applicability of the vesicles. Here, the dynamics of lipid molecules within vesicles of approximately 10 nm in diameter are monitored using coarse-grained molecular dynamics (CGMD) simulations (MARTINI). The vesicles consist of (i) 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), as well as 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), each in a mixture with cholesterol (CHOL) (POPC/CHOL, DOPC/CHOL), (ii) POPC/DOPC mixture and (iii) pure POPC and DOPC. Vesicles are generated by placing the individual lipids spatially separated by an area-per-lipid based distance onto the sphere using an arc correction, implemented into ArcVes, which is an in-house developed software. This protocol leads to the generation of vesicles which remain stable during extended MD simulations. In the presence of cholesterol, the mixed vesicles show time-dependent changes such as a decrease in both radii and membrane thickness. The number of lipids in each leaflet of the mixed vesicles reveals a trend of POPC and DOPC flipping to the outer leaflet and of cholesterol to the inner leaflet. This leads to a cholesterol rich path-formation in both the inner and outer leaflets. The diffusivity of the lipids in the vesicles mixed with cholesterol is lower than the diffusivity of the pure vesicles, similar to the observation for the respective flat membrane patches. In general, the diffusivity of the lipids is larger in the outer leaflets than the inner leaflets.
Citation: Chia-Wen Wang, Meng-Han Lin, Wolfgang B. Fischer. Cholesterol affected dynamics of lipids in tailor-made vesicles by ArcVes software during multi micro second coarse grained molecular dynamics simulations[J]. AIMS Biophysics, 2023, 10(4): 482-502. doi: 10.3934/biophy.2023027
Related Papers:
[1]
Yong Ding, Jian-Hong Liu .
The signature lncRNAs associated with the lung adenocarcinoma patients prognosis. Mathematical Biosciences and Engineering, 2020, 17(2): 1593-1603.
doi: 10.3934/mbe.2020083
[2]
Rui Huang, Xiwen Liao, Qiaochuan Li .
Integrative genomic analysis of a novel small nucleolar RNAs prognostic signature in patients with acute myelocytic leukemia. Mathematical Biosciences and Engineering, 2022, 19(3): 2424-2452.
doi: 10.3934/mbe.2022112
[3]
Siqi Hu, Fang Wang, Junjun Yang, Xingxiang Xu .
Elevated ADAR expression is significantly linked to shorter overall survival and immune infiltration in patients with lung adenocarcinoma. Mathematical Biosciences and Engineering, 2023, 20(10): 18063-18082.
doi: 10.3934/mbe.2023802
[4]
Wei Niu, Lianping Jiang .
A seven-gene prognostic model related to immune checkpoint PD-1 revealing overall survival in patients with lung adenocarcinoma. Mathematical Biosciences and Engineering, 2021, 18(5): 6136-6154.
doi: 10.3934/mbe.2021307
[5]
Jinqi He, Wenjing Zhang, Faxiang Li, Yan Yu .
Development of metastasis-associated seven gene signature for predicting lung adenocarcinoma prognosis using single-cell RNA sequencing data. Mathematical Biosciences and Engineering, 2021, 18(5): 5959-5977.
doi: 10.3934/mbe.2021298
[6]
Yuan Yang, Lingshan Zhou, Xi Gou, Guozhi Wu, Ya Zheng, Min Liu, Zhaofeng Chen, Yuping Wang, Rui Ji, Qinghong Guo, Yongning Zhou .
Comprehensive analysis to identify DNA damage response-related lncRNA pairs as a prognostic and therapeutic biomarker in gastric cancer. Mathematical Biosciences and Engineering, 2022, 19(1): 595-611.
doi: 10.3934/mbe.2022026
Yi Liu, Long Cheng, Xiangyang Song, Chao Li, Jiantao Zhang, Lei Wang .
A TP53-associated immune prognostic signature for the prediction of the overall survival and therapeutic responses in pancreatic cancer. Mathematical Biosciences and Engineering, 2022, 19(1): 191-208.
doi: 10.3934/mbe.2022010
[9]
Ji-Ming Wu, Wang-Ren Qiu, Zi Liu, Zhao-Chun Xu, Shou-Hua Zhang .
Integrative approach for classifying male tumors based on DNA methylation 450K data. Mathematical Biosciences and Engineering, 2023, 20(11): 19133-19151.
doi: 10.3934/mbe.2023845
[10]
Pingping Song, Jing Chen, Xu Zhang, Xiaofeng Yin .
Construction of competitive endogenous RNA network related to circular RNA and prognostic nomogram model in lung adenocarcinoma. Mathematical Biosciences and Engineering, 2021, 18(6): 9806-9821.
doi: 10.3934/mbe.2021481
Abstract
In biotechnology, as well as in nanomedicine, vesicular structures are essential features of the cellular life cycle. Lipid compositions define the applicability of the vesicles. Here, the dynamics of lipid molecules within vesicles of approximately 10 nm in diameter are monitored using coarse-grained molecular dynamics (CGMD) simulations (MARTINI). The vesicles consist of (i) 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), as well as 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), each in a mixture with cholesterol (CHOL) (POPC/CHOL, DOPC/CHOL), (ii) POPC/DOPC mixture and (iii) pure POPC and DOPC. Vesicles are generated by placing the individual lipids spatially separated by an area-per-lipid based distance onto the sphere using an arc correction, implemented into ArcVes, which is an in-house developed software. This protocol leads to the generation of vesicles which remain stable during extended MD simulations. In the presence of cholesterol, the mixed vesicles show time-dependent changes such as a decrease in both radii and membrane thickness. The number of lipids in each leaflet of the mixed vesicles reveals a trend of POPC and DOPC flipping to the outer leaflet and of cholesterol to the inner leaflet. This leads to a cholesterol rich path-formation in both the inner and outer leaflets. The diffusivity of the lipids in the vesicles mixed with cholesterol is lower than the diffusivity of the pure vesicles, similar to the observation for the respective flat membrane patches. In general, the diffusivity of the lipids is larger in the outer leaflets than the inner leaflets.
1.
Introduction
The susceptible-infective-removed-susceptible (SIRS) model is one of the most popular epidemic models, in which total host population is divide into three classes called susceptible(S), infective(I) and removed(R). Since the Kermack and McKendrick pioneering work [1], a large number of studies have been carried out on the model of continuous infectious diseases attempting to gain a better understanding of diseases transmission, especially for the control policies and dynamics [2,3,4,5]. Recent years, many researchers have studied the epidemic with SIRS models such as COVID-19(2019), hand-foot-and-mouth diseases(HFMD) and other infectious diseases recently. Although scholars have flocked to study many properties of the SIRS models, there are still many aspects that deserve further study, such as the threshold dynamics of SIRS model with vaccination and optimal control using dynamic programming method.
In studying the dynamics of the disease, one of the most important concepts is the basic reproduction number R0. It is epidemiologically defined as the expected number of secondary cases produced by a typical infectious individual during its entire infectious period in a completely susceptible host population. To investigate the global behavior of the prevalence of infectious diseases, the stability analysis of equilibrium for epidemic models have been carried out (see [6,7] and the references therein). For ordinary differential SIRS models, Mena-Lorca and Hethcote considered several kinds of SIRS epidemic models and a threshold parameters were also found in [6] to determine whether the disease dies out or approaches to an endemic equilibrium, we also note that the global stability of SIRS epidemic model have already been investigated and the threshold theorems are well obtained. In the subsequent studies, it was generally believed that the introduction of age-structure into epidemic models is necessary and reasonable. Since the fact that the age structure of a population affects the dynamics of disease transmission was recognized, various age-structured epidemic models have been investigated by many authors. Some theoretics on the transmission dynamics of age-structured epidemic models have been developed in [8,9,10,11]. He et al. [12] studied the optimal birth control of age-dependent competitive species. Yang et al. [13] proposed epidemic model with age-since-infection and diffusion, the next generation operator R is also given as the basic regeneration numbers. In paper [14], Yang et al. also studied the threshold dynamics by the method of Lyapunov functionals for an age-of-infection epidemiological model with nonlinear incidence rate. Chekroun et al. [15] showed the stability of the equilibrium by using Lyapunov functions. Lan et al. [16] just studied the impact of hospital resources and environmental perturbation to the dynamics of SIRS model. But the effect of vaccination on threshold dynamics was not considered in these articles.
In fact, we notice that vaccination has a significant impact on the thresholds dynamics of many infectious diseases, such as COVID-19, HPV, etc. Obviously, vaccination of susceptible individuals can greatly reduce the probability that susceptible individuals become infected. Thus, decision makers hope to achieve the optimal proportion of vaccines with the lowest cost. Optimal control is a promising strategy to control disease outbreaks. The main purpose of using the control in some diseases is to search for the most effective one that reduces the infection to a minimum level while minimizing the cost of applying control measures. In papers [7,9,17], the strategies such as health promotion campaigns or lockdown policies can contribute to reducing the disease transmission. In paper [18], Mu and Zhang investigated the near-optimal control of SIRS multi-strain epidemic model with age-structure using Pontryagin maximum principle. Bolzoni et al. [19] investigated optimal control of epidemic size and duration with limited resources. Zhou et al. [20] researched optimal isolation strategies of emerging infectious diseases with limited resources. In papers [21,22,23], some measures correspond to screening and quarantining of infected were implemented to achieve the purpose of controlling the disease.
However, we note that the most of these articles by using maximum principle to show optimal controls and only present necessary conditions. Defectively, the method of maximum principle only addresses the optimal control problem with initial value case of t0=0. In fact, vaccination can be administered at any time t=t0 in the control of the disease. Surpassingly, the dynamic programming method can give the necessary and sufficient conditions for the existence of optimal control via the HJB equation at any initial times t=t0 and initial states. In other words, the dynamic programming method is an extension of the maximum principle, namely, the dynamic programming method is equivalent to maximum principle method when t0=0.
Overall, the threshold dynamics results about age-structured SIRS epidemic models with the proportion of the vaccine injections are comparatively scarce, as well as few articles show that dynamic programming approach to solve the optimal control problem. Thus, in this paper, taking into account for the effect of the proportion of the vaccine injections, we incorporate the vaccine injection rate in solving the threshold of age-structured SIRS epidemic models. And vaccination optimal control is also obtained through dynamic programming approach.
The layout of this paper is as follows: in section 2, an age-structure SIRS epidemic model is reproduced and new age-structure SIRS epidemic model with vaccination is given; in section 3, we give the basic regeneration number and discuss the threshold dynamics; in section 4, we formulate and solve the corresponding optimal control problem. The existence of optimal control is proved by HJB equation and the expression of optimal control is present. in section 5, some numerical simulations are performed to demonstrate the theoretical results. Brief conclusions are given in section 6.
2.
Model and preliminaries
We reproduce the age-structured SIRS epidemic model based on paper [24] as follows
where k1,k2∈[0,1] present the weight of the infected and recovered, respectively. All the parameters are related to age a, positive and bounded, These parameters and their meaning are listed in the Table 1, where η(a)=σ(a)+α(a), and δ(a)=μ(a)+ϕ(a)+σ(a)+α(a).
Table 1.
Parameter values used in numerical simulations.
Notation
Biological meanings
S(a,t)
Densities of susceptible individuals of age a at time t
I(a,t)
Densities of infective individuals of age a at time t
R(a,t)
Densities of recovery individuals of age a at time t
As we all know that the age density distribution function of the total population reaches a stable state, i.e
S(a,t)+I(a,t)+R(a,t)=P(a).
(2.3)
We focus on threshold dynamics on an age-structured SIRS epidemic model with vaccination. The spread of many epidemics are age-related, such as Human Papilloma Virus (HPV) spread across the population after ten years olds usually. And the World Health Organization believes that the most appropriate age for HPV vaccination was 11-12 years old. Therefore, it makes sense to study the vaccination is only given for a certain age to age-structured SIRS epidemic model.
The introduced vaccination rate u(A0) satisfies 0≤u(A0)≤1 and is only given for a certain age and susceptible individuals of age A0 is inoculated at some point. The disease does not spread among people younger than age A0, that is, there are no patients before the vaccination age A0. Therefore, we can get that the system satisfied by the susceptible, infective, and recovering individuals based on model (2.1) as following
where S(A−0,t) and I(A−0,t) are the densities of susceptible and recovering individuals in the left neighborhood of age A0, k(a) is a positive continuous function. The vaccination rate u(A0) is a constant and k(a) is a nonnegative continuous function defined on 0≤a≤A, where A is the maximum age.
In this paper, unless otherwise specified. The H represents the Hilbert space, with inner product
⟨φ1,ϕ1⟩=∫A0φ1(a)ϕ1(a)da,φ1,ϕ1∈H,
(2.5)
and we consider the state equations on space L1(0,A;H3), when φ=(φ1,φ2,φ3) and ϕ=(ϕ1,ϕ2,ϕ3), φi,ϕi∈H, the definition of inner product on L1(0,A;H3) as
In order to obtain the endemic disease equilibrium solution, S(a,t)=S(a),I(a,t)=I(a) and R(a,t)=R(a) are substituted into the model (2.4) obtained
R(a)=0,0<a<A0,R(a)=u(A0)P(a),A0<a<A.
(3.2)
And we know when R0(a)=u(A0)P(a), the R(a,t)=u(A0)P(a), and when the population reaches a dynamic equilibrium P(a), Combining the equations (2.3) and (3.2), and substitute into (2.4) for the equation that satisfies I(a,t). From (2.3), it shows that
It is quite obvious the right-hand side of the equation (3.8) is the minus function about D∞, thus, we can define
R0=(1−u(A0))∫A0∫aA0k(τ)p∞(τ)⋅exp(−∫aτδ(θ)dθ)dτda.
(3.9)
When R0<1, (3.7) has only a unique zero solution and no positive root, when R0>1, the (3.7) formula has zero root and the only positive root I∗>0 satisfies (3.8) formula, thus there are only positive root E∗=(S∗,I∗,R∗) of (2.4). For ease of analysis, we converse i(a,t)=I(a,t)(1−u(A0))P, equation (2.4) about I(a,t) can be transformed into
From the transformation, we can see that the existence and stability of the positive equilibrium of model (2.4) is equivalent to positive equilibrium of (3.10) model. R0 can be used as a threshold for the existence or absence of endemic diseases. It is also a threshold for the disappearance of diseases. Let ℏ=1−1R0 and we call ℏ vaccination elimination point.
Similar to the [25] method, the following conclusion can be obtained. Before discussing the global stability of equilibriums, we first propose the lemma 3.1, and the similar proof process have been presented in [26].
Lemma 3.1.Let i1(a,t) and i2(a,t) are the solutions of the (3.10) equations, satisfying i1(a,0)=i10(a)≤i2(a,0)=i20(a), where i10(a),i20(a)∈K, then there are following properties
(1) ij(a,t)∈K,j=1,2;
(2) i1(a,t)≤i2(a,t),t≥0;
(3) When 0<ξ<1, the solution i(a,t) of (3.10) that the initial value satisfies the i(a,0)=i10(a)ξ, then it must satisfy the i1(a,t)ξ≤i(a,t),
where the K is denoted
K={f(a,t)|f(a,t)∈C,f(0,t)=0,0≤f(a,t)≤1}.
Theorem 3.2.When R0<1, the model (3.10) has only a disease-free equilibrium, which is globally asymptotically stable, while when R0>1, the model (3.10) also has a unique endemic equilibrium E∗, which is globally asymptotically stable and the disease-free equilibrium is unstable.
Proof. We only consider the asymptotic behavior of solutions where the initial values of the model (3.10) satisfy 0≤i(a,0)≤1 in (A0,A). The theorem of the global stability is given for solutions.
It is not difficult to see if the model (3.10) has a positive equilibrium equivalent to the (3.5) has a positive root. Since the right of the (3.5) is reduced about D(t) and tends to zero when the D(t) tends to positive infinity. There is no positive real root to the (3.5) when R0<1, and there is a unique positive real root when R0>1, i.e. the model (2.4) has only a unique positive equilibrium when R0>1 and, instead, no positive equilibrium when R0<1.
When R0<1, in order to study the global stability of the disease-free equilibrium, expression of the solution of (3.10) about I(a,t) is obtained by using the characteristic line method
So we obtain χ(t)≤R0max{χn−1,χn} when nA≤t≤(n+1)A. Because of R0<1, we also obtain χn<Ln=χn−1, i.e. χn<χn−1, χn is monotonous and reduced and χn<χn−1R0, from this delivery limn→+∞χn=0. Considering above conditions and (3.13) can also derive limt→+∞i(a,t)=0. Finally, it is concluded from (3.13) that the disease-free equilibrium of the model (3.10) is stable R0<1.
When R0>1, substituting i(a,t)=i(a) into the (3.10), the only endemic equilibrium of (3.10) is given by
When A0≤a<A, the ∀i(a,t) satisfy i(a,0)>0, and χ(t)>0 can be derived from (3.11) and (3.14) when 0≤t≤A, and we also know that i(a,A)>0(A0<a<A) by (3.13). Repeat this process then the χ(t)>0 is obtained when 0≤t≤3A, let χ0=minA≤t≤2Aχ(t), χ∗=maxA≤t≤2Aχ(t), the following is got as
Comparing the forms of (3.16) and (3.17), we can find a positive number 0<ξ<1 such that i(a,2A)≥ξi(a), namely, the following formula holds for the solution i(a,t) of (2.4)
ξ⋅i(a)≤i(a,2A)≤1.
The Lemma 3.1 shows that the solutions iξ(a,t) and i1(a,t) with initial values of ξ⋅i(a,0) and 1 respectively when t=0, must have
iξ(a,t)≤i(a,t+2A)≤i1(a,t).
Based on the fact that i(a) is the equilibrium of (3.10), then the ξi(a)≤iξ(a,t), i.e. iξ(a,t) is a monotone increasing function about t. Since i1(a,t) is a monotone subtraction function about t, and the corresponding χξ(t) and χ1(t) are monotone increasing and monotone decreasing, respectively, and there are limt→∞χξ(t)=χ0ξ≤limt→∞χ1(t)=χ01. From (3.14), the equation has only unique normal number solutions when R0>1, i.e. χ0ξ=χ01=ξ∗, the χ(t) corresponding to the i(a,t) also obtain limt→∞χ(t)=ξ∗. Reuse (3.13) when t>A0, can be obtained limt→∞i(a,t)=i(a), i.e. the local disease equilibrium of (3.10) is globally attractive. The stability of the local disease equilibrium can be obtained by linearization. Therefore, the local disease equilibrium is globally asymptotically stable.
Remark 3.3.3.2 shows that the epidemic disease will die out if R0<1, while the disease will be prevalent if R0>1. This implies that R0 is the threshold of model (2.4).
4.
Optimal control strategies
4.1. Formulation of the optimal control problem
In section 3, we consider the effect of a age determined vaccine immunization of the system on threshold dynamics, which is one of the cases in which vaccines are added to all susceptible persons in the system. In this section, we focus on the optimal control of (4.1) on the basis of age-structured SIRS by implementing proportion of vaccination u(a,t) as control strategy at all ages of system, and vaccination rate u(a,t) satisfies 0≤u(a,t)≤1. We aim to minimize the infected population, while keeping the cost of applying such control strategy at the minimum level.
Let ((a,t0),x(a,t0))∈[0,T)×(H×H×H), and consider the following control system over on (a,t)∈(0,A)×[t0,T]:
where S(a,t0)=St0(a),I(a,t0)=It0(a),R(a,t0)=Rt0(a) is varying initial time value and denote the initial state x(a,t0) as xt0(a). Here the control {u(⋅,⋅):(0,A)×[t0,T]→U}, and u(⋅,⋅) is measurable. We obtain the schematic diagram of model (4.1) with introducing vaccination control (see Figure 1).
The object is to design the optimal controller u∗(a,t),(a,t)∈(0,A)×[t0,T] for the system (4.1) which minimizes the following performance index during the finite time interval [t0,T][27,28,29,30]. Our goal is to minimize the total number of the infected and susceptible individuals by using minimal control efforts. The cost function as follows
where I(a,t):(a,t)∈[0,A]×[t0,T] and we denote x=(S,I,R),xt0∈(H×H×H), the terminal value is V(x,T)=∫A0h(I(a,T))da. Let L(t,x,u):=τ1I(a,t)+τ2u(a,t)+12τ3u2(a,t), where τ1, τ2 and τ3 are the positive weighting factors, representing the cost per unit time of the components I(a,t), u(a,t) and u2(a,t), respectively. In particular, ∫Tt0∫A0τ1I(a,t)dadt is the cost that infected individual creates for the society due to standard medical care, not including the vaccination treatment u(a,t). ∫Tt0∫A0(τ2u(a,t)+12τ3u2(a,t))dadt is the cost of treating infected individuals. h:H×H×H→H are measurable and h(I(a,T)) denotes the total number disease in population depend on the age at the time T. The value function as
The function V(⋅,⋅) will play an important role in obtaining optimal control, thus, we would like to study V(⋅,⋅) in great detail. Let (x∗,u∗) is optimal pair, where x∗=(S∗,I∗,R∗). Noting that the initial time (t=0) and the initial state (x(a,0)=x0) are fixed in the model [18]. The h(I(a,T)) satisfies
∫A0∣h(I(a,T))−h(˜I(a,T))∣da≤C‖I(a,T)−˜I(a,T)‖C,
(4.4)
where C is a positive constant, ‖⋅‖ represents norm on space H and I(a,T), ˜I(a,T)∈H. The following hypothesis are put forward and need to be satisfied.
D1(U,d) is a separable metric space and T>0.
D2 The control set U is bounded convex
D3 Different controls correspond to same terminal value.
Then we study dynamic programming method, to give optimal strategy of model (4.1). The basic idea of this approach applied to optimal controls is to consider a family of optimal control problems with different initial times (t=t0≥0) and initial states (x(a,t0)=xt0(a)) [see [31,32,33]], to establish relationships among these problems through the HJB equation. Firstly, we present the Bellman's principle of optimality as follows.
Theorem 4.1.Let (D1)-(D3) hold. Then for any (t0,xt0(a))∈[0,T)×H×H×H,
The equation (4.5) is dynamic programming equation, which gives a relationship among globally and locally optimal via value function.
We would like to further study the (4.5), trying to obtain an equation for V(t0,xt0(a)) simpler form. The following results give a partial differential equation that a continuously differentiable value function ought to satisfy. We know the V(t0,xt0(a)) is a continuously differentiable value function via 4.1 and basic theorem of calculus.
Proof. Fix (t0,xt0(a))∈[0,T)×H×H×H and u∈U. Let x(⋅,⋅) be the state trajectory corresponding to the control u(⋅,⋅)∈Uω[s,T] with u(a,t)≡u. On the one hand, we can write according to (4.5)
Combining (4.14) and (4.17), we obtain our results.
Definition 4.3.[34] A function v∈C([0,T]×H×H×H) is called a viscositysolution of (4.9) if
v(T,x(a,T))≤h(x(a,T)),∀x(a,T)∈H×H×H,
(4.18)
and for any ψ∈C1([0,T]×R), whenever v−ψ attains a local maximum at (t,x)∈[0,T)×H×H×H, we have
−ψt(t,x)+supu∈UH(t,x,u,−ψx(t,x))≤0.
(4.19)
The function v∈C([0,T]×H×H×H) is called a viscositysolution of (4.9) if in (4.18)−(4.19) the inequalities ≤ are changed to ≥ and localmaximum is changed to localminimum. In the case that v is both a viscosity subsolution and supersolution of (4.9), it is called a viscosity solution of (4.9).
Theorem 4.4.Let (D1)-(D3) hold. Then the value function of V(t,x) satisfies
Thus, V(t,x(a,t)) satisfies (4.21). We can know that V(t,x(a,t)) is unique viscosity of the HJB equation via the same method as [[34], Theorem 2.5].
4.2. Optimal control
In this subsection, we study the existence of the optimal control pair in finite time for the system (4.1) and construct the Hamiltonian of the optimal control problem.
Theorem 4.5.There exists an optimal control u∗∈U and a corresponding optimal state S∗,I∗,R∗ such that
Proof. To show the existence of the optimal control for the problem under consideration, we use the result in [35,36]. We notice that the state and control variables are non-negative, and the control set U, by definition, is closed and bounded. The optimal system is bounded, which ensures the compactness needed for the existence of an optimal control. Further, the integrand ∫Tt0∫A0(τ1I(a,t)+τ2u(a,t)+12τ3u2(a,t))dadt is convex on the control set U due to the biquadratic and quadratic nature of control variable u(a,t), which represents vaccination proportion. In addition, there exists a constant ν>1 and positive numbers κ1, κ2 such that
τ1I(a,t)+τ2u(a,t)+12τ3u2(a,t)≥κ1|u(a,t)|ν−κ2.
We complete the existence of the optimal control of state variables.
Theorem 4.6.Let u∗(a,t) be optimal control variable, S∗(a,t), I∗(a,t) and R∗(a,t) be corresponding optimal state variable of model (4.1). The corresponding optimal control is given as follows:
Proof. We define the Hamiltonian function H(t,x,u(a,t),Vx),t∈[0,T],u(a,t)∈Uad,x,Vx(x,t)∈H×H×H and by constructing the Hamiltonian function and HJB equation to characterize this optimal control u(a,t) as follows:
Let u∗(a,t) be given optimal control and S∗(a,t), I∗(a,t) and R∗(a,t) be corresponding optimal state variable of model (4.1). By the Hamiltonian function (4.22), optimal control u(a,t)=u∗(a,t) is obtained
u∗(a,t)=−1τ3S(a,t)∗(VR(x,t)−VS(x,t))+τ.
(4.23)
From the properties of the control set with the findings, then there is exists optimal control
Remark 4.7.If we could obtain the value function V(t,x(a,t)) by solving the HJB equation, then we could, at least formally, construct an optimal pair for each of optimal control problem for different initial values. The principle involves the following steps to solve optimal control problem
Step.1. Solve the HJB equation (4.10) to find the V(T,x(a,t));
Step.2. Find u∗(a,t) through Hamiltonian function (4.22);
Step.3. Combine the u∗(a,t) to solve model (4.1) to get the optimal pair (x∗(a,t),u∗(a,t)).
5.
Numerical simulation
This section aims to illustrate the effectiveness of our theoretical results obtained in previous sections.
5.1. Numerical simulation of the stability of equilibria
In the numerical simulation, the parameters are given in Table 2. We assume that (a,t)∈[0,80]×[0,500],Δt=0.1,Δa=1, take p∞=1.0,ϱ=0.2,δ=0.1827,k=0.2, then R0≈2.69>1. The parameters that we partially need to assume are reasonable according to the study of Yang et al. [14] has shown that the range of the basic reproduction number of the pandemic influenza in 1.0−4.0. We see from 3.2 the E∗ is globally asymptotically stable. Actually, as showing in the following Figure 2 (a), the density of the infected individual I(a,t) tends to be a positive constant over time.
Table 2.
Parameter values used in numerical simulations.
Take ϱ=0.1,δ=0.1827,k=0.1, then R0≈0.78<1, thus we can know the E0 is globally asymptotically stable, as showing in the following Figure 2 (b), the density of the infected individual I(a,t) tends to zero over time. These are the same conclusions as 3.2.
We choose δ=0.1827, 0.2027 and 0.2227 to study the effect on the infected population and all the other parameters are fixed as in Table 2. We calculate the corresponding basic reproduction number R0≈0.9427, 0.8979 and 0.7802 when R0<1, respectively. In Figure 3 (a), we can see that the extinction time of decreases with the increase of infected removal rate δ. In order to simulate the cases of R0>1, we calculate the corresponding basic reproduction number R0≈1.3468, 2.2441 and 2.6930, respectively. The simulation result is presented in Figure 3 (b). Synthetically, it is obvious to see from Figure 3 that the peak of I(a,t) decreases significantly as infected removal rate δ of age a goes up.
Figure 3.
(a): the path of I(a,t) under different δ when R0<1 at age a=30; (b): the path of I(a,t) under different δ when R0>1 at age a=30.
In this section, we show the result of optimal control using a numerical example. The parameters are given in Table 2. (a,t)∈[0,80]×[0,500],Δt=0.1,Δa=1 by performing some numerical simulations. We take τ1=11.7, τ2=3.4 and τ3=12. The corresponding paths of control and infective populations I(a,t) are plotted in Figure 4, Figure 5 (a), Figure 6 and Figure 7 (a) show the section view at age a=20 and a=30 of the image variation tendency of I(a,t) with and without control and the tendency of corresponding optimal control u(a,t). We observe that the comprehensive use of control u(a,t) works better than no control applied during the course of disease. The effect of applying control reduces the number of infective populations I(a,t) dramatically. In short, control variable u(a,t) play an important role in the control of the disease. At the same time, we can also observe a reduced epidemic time of population infection.
Figure 4.
When R0>1, (a): the paths of I(a,t) with and without control u(a,t) at age a=20; (b): The paths of I(a,t) with and without control u(a,t) at age a=30.
Figure 6.
When R0<1, (a): the paths of I(a,t) with and without control u(a,t) at age a=20; (b): The paths of I(a,t) with and without control u(a,t) at age a=30.
We note that a significant decline in the number of infected persons I(a,t) during vaccine control of susceptible at age a=20 and a=30 from these following figures.
When R0>1, I(a,t) trends is regard to the optimal controls, where the upper-level image is the evolution of I(a,t) without control and the lower-level image is the evolution of I(a,t) with control on in Figure 4. In Figure 5 (a) shows the optimal control pathes at ages a=20 and a=30, repectively. Eventually the disease tends to extinction with control. Compared with the path of non-vaccinated I(a,t), the epidemic time of I(a,t) with vaccination is significantly shortened, as well as the density of infected population verge to zero eventually.
When R0<1, I(a,t) trends is regard to the optimal controls, where the upper-level image is the evolution of I(a,t) without control and the lower-level image is the evolution of I(a,t) with control in Figure 6. In Figure 7 (a) shows the optimal control pathes at ages a=20 and a=30, repectively. Eventually the disease tends to extinction with control. Compared with Compared with the path of non-vaccinated I(a,t), the epidemic time of I(a,t) with vaccination is significantly shortened. Combined with Figure 4 and the corresponding conclusions, it is clear that vaccination have a certain degree of hindrance to the transmission of epidemics. We also observe the Figure 5 (b) and Figure 7 (b) that vaccine control varies for different ages, mainly adolescents and middle age, and this control strategy is consistent with the actual control of many epidemics such as pneumonia, influenza and so on. The pneumococcal vaccine, COVID-19 vaccine and influenza vaccine are not suitable for the elderly. Overall, population age has a significant effect on the implementation of vaccine control strategies for many common diseases.
6.
Concluding remarks
In this paper, we recurred an age-structured SIRS epidemic model (2.1) with vaccination control, and obtained the condition for disease extinction and persistence. The research results showed that when R0<0, the disease extincted; when R0>0, the disease persisted. The conditions of disease extinction and permanence were given in 3.2. And we verified the resluts by numerical simulation. However, since the bifurcation phenomenon will occur when R0=0, this part will also be an important point for our further study of the age-structured SIRS epidemic system dynamics. Besides, we also studied the control problem of SIRS epidemic model with age structure. By proving the existence of viscosity solution, we obtain the existence of optimal control. Then, we obtained the implicit expression for optimal control and HJB equation. Finally, by using Milstein method, the optimal proportion of vaccination is obtained. In addition, we plan that in the next step, traveling wave solutions of an age-structured SIRS epidemic model will become one of our main research contents. Wu et al. [37] studied the existence and non-existence of traveling wave solutions for a diffusive age-structured SIR epidemic model. Their work provides some insights on how to deal with the high dimensional age-structured epidemic models.
Acknowledgments
The authors thank the editor and referees for their careful reading and valuable comments. The research was supported in part by the Ningxia Natural Science Foundation (No. 2020AAC03067).
Conflict of interest
The authors have no conflict of interest.
Acknowledgments
WBF thanks the Ministry of Science and Technology Taiwan (MOST-107-2112-M-010-002-MY3 and MOST-110-2112-M-A49A501) for financial support.
van Meer G, Voelker DR, Feigenson GW (2008) Membrane lipids: where they are and how they behave. Nat Rev Mol Cell Biol 9: 112-124. https://doi.org/10.1038/nrm2330
[3]
Akbarzadeh A, Rezaei-Sadabady R, Davaran S, et al. (2013) Liposome: classification, preparation, and applications. Nanoscale Res Lett 8: 102. https://doi.org/10.1186/1556-276X-8-102
Lawrence RE, Zoncu R (2019) The lysosome as a cellular centre for signalling, metabolism and quality control. Nat Cell Biol 21: 133-142. https://doi.org/10.1038/s41556-018-0244-7
[6]
Miorin L, Romero-Brey I, Maiuri P, et al. (2013) Three-dimensional architecture of thick-borne encephalitis virus replication sites and trafficking of the replication RNA. J Virol 87: 6469-6481. https://doi.org/10.1128/JVI.03456-12
[7]
Rideau E, Dimova R, Schwille P, et al. (2018) Liposomes and polymersomes: a comparative review towards cell mimicking. Chem Soc Rev 47: 8572-8610. https://doi.org/10.1039/c8cs00162f
Sercombe L, Veerati T, Moheimani F, et al. (2015) Advances and challenges of liposome assisted drug delivery. Front Pharmacol 6: 286. https://doi.org/10.3389/fphar.2015.00286
[10]
Bulbake U, Doppalapudi S, Kommineni N, et al. (2017) Liposomal formulations in clinical use: an updated review. Pharmaceutics 9: E12. https://doi.org/10.3390/pharmaceutics9020012
Gudlur S, Sandén C, Matoušková P, et al. (2015) Liposomes as nanoreactors for the photochemical synthesis of gold nanoparticles. J Colloid Interface Sci 456: 206-209. https://doi.org/10.1016/j.jcis.2015.06.033
[13]
Marrink SJ, Mark AE (2003) Molecular dynamics simulation of the formation, structure, and dynamics of small phospholipid vesicles. J Am Chem Soc 125: 15233-15242. https://doi.org/10.1021/ja0352092
[14]
Reynwar BJ, Illya G, Harmandaris VA, et al. (2007) Aggregation and vesiculation of membrane proteins by curvature-mediated interactions. Nature 447: 461-464. https://doi.org/10.1038/nature05840
[15]
Bryer AJ, Reddy T, Lyman E, et al. (2022) Full scale structural, mechanical and dynamical properties of HIV-1 liposomes. PLoS Comput Biol 18: e1009781. https://doi.org/10.1371/journal.pcbi.1009781
[16]
Nawae W, Hannongbua S, Ruengjitchatchawalya M (2014) Defining the membrane disruption mechanism of kalata B1 via coarse-grained molecular dynamics simulations. Sci Rep 4: 3933. https://doi.org/10.1038/srep03933
[17]
Kawamoto S, Klein ML, Shinoda W (2015) Coarse-grained molecular dynamics study of membrane fusion: Curvature effects on free energy barriers along the stalk mechanism. J Chem Phys 143: 243112. https://doi.org/10.1063/1.4933087
Parton DL, Tek A, Baaden M, et al. (2013) Formation of raft-like assemblies within clusters of influenza hemagglutinin observed by MD simulations. PLoS Comput Biol 9: e1003034. https://doi.org/10.1371/journal.pcbi.1003034
[20]
Louhivuori M, Risselada J, van der Giessen E, et al. (2010) Release of content through mechano-sensitive gates in pressurized liposomes. Proceed Natl Acad Sci USA 107: 19856-19860. https://doi.org/10.1073/pnas.1001316107
[21]
Dieluweit S, Csiszár A, Rubner W, et al. (2010) Mechanical properties of bare and protein-coated giant unilamellar phospholipid vesicles. A comparative study of micropipet aspiration and atomic force microscopy. Langmuir 26: 11041-11049. https://doi.org/10.1021/la1005242
Dao TPT, Fauquignon M, Fernandes F, et al. (2017) Membrane properties of giant polymer and lipid vesicles obtained by electroformation and pva gel-assisted hydration methods. Colloid Surface A 533: 347-353. https://doi.org/10.1016/j.colsurfa.2017.09.005
[25]
Steinkühler J, De Tillieux P, Knorr RL, et al. (2018) Charged giant unilamellar vesicles prepared by electroformation exhibit nanotubes and transbilayer lipid asymmetry. Sci Rep 8: 11838. https://doi.org/10.1038/s41598-018-30286-z
[26]
Roy B, Guha P, Bhattarai R, et al. (2016) Influence of lipid composition, pH, and temperature on physicochemical properties of liposomes with curcumin as model drug. J Oleo Sci 65: 399-411. https://doi.org/10.5650/jos.ess15229
[27]
de Vries AH, Mark AE, Marrink SJ (2004) Molecular dynamics simulation of the spontaneous formation of a small DPPC vesicle in water in atomistic detail. J Am Chem Soc 126: 4488-4489. https://doi.org/10.1021/ja0398417
[28]
Shinoda W, DeVane R, Klein ML (2010) Zwitterionic lipid assemblies: molecular dynamics studies of monolayers, bilayers, and vesicles using a new coarse grain force field. J Phys Chem B 114: 6836-6849. https://doi.org/10.1021/jp9107206
[29]
Shillcock JC (2012) Spontaneous vesicle self-assembly: a mesoscopic view of membrane dynamics. Langmuir 28: 541-547. https://doi.org/10.1021/la2033803
[30]
Qi Y, Ingólfsson HI, Cheng X, et al. (2015) CHARMM-GUI Martini maker for coarse-grained simulations with the martini force field. J Chem Theory Comput 11: 4486-4494. https://doi.org/10.1021/acs.jctc.5b00513
[31]
Hsu PC, Bruininks BMH, Jefferies D, et al. (2017) CHARMM-GUI Martini maker for modeling and simulation of complex bacterial membranes with lipopolysaccharides. J Comput Chem 38: 2354-2363. https://doi.org/10.1002/jcc.24895
[32]
Damre M, Marchetto A, Giorgetti A (2019) MERMAID: dedicated web server to prepare and run coarse-grained membrane protein dynamics. Nucleic Acids Res 47: W456-W461. https://doi.org/10.1093/nar/gkz416
[33]
Risselada HJ, Mark AE, Marrink SJ (2008) Application of mean field boundary potentials in simulations of lipid vesicles. J Phys Chem B 112: 7438-7447. https://doi.org/10.1021/jp0758519
[34]
Sharma SD, Kim BN, Stansfeld PJ, et al. (2015) A coarse grained model for a lipid membrane with physiological composition and leaflet asymmetry. PLoS One 10: e0144814. https://doi.org/10.1371/journal.pone.0144814
[35]
Giuliari B, Kösters M, Zhou J, et al. The Vesicle Builder–A Membrane Packing Algorithm for the CELLmicrocosmos MembraneEditor (2020). https://doi.org/10.2312/molva.20201096
[36]
Jo S, Kim T, Iyer V, et al. (2008) CHARMM-GUI: a web-based graphical user interface for CHARMM. J Comput Chem 29: 1859-1865. https://doi.org/10.1002/jcc.20945
[37]
Brooks BR, Brooks III CL, Mackerell Jr. AD, et al. (2009) CHARMM: the biomolecular simulation program. J Comput Chem 30: 1545-1614. https://doi.org/10.1002/jcc.21287
[38]
Lee J, Cheng X, Swails JM, et al. (2016) CHARMM-GUI input generator for NAMD, GROMACS, AMBER, OpenMM, and CHARMM/OpenMM simulations using the CHARMM36 additive force field. J Chem Theory Comput 12: 405-413. https://doi.org/10.1021/acs.jctc.5b00935
[39]
Rivel T, Ramseyer C, Yesylevskyy S (2019) The asymmetry of plasma membranes and their cholesterol content influence the uptake of cisplatin. Sci Rep 9: 5627. https://doi.org/10.1038/s41598-019-41903-w
[40]
Wilson KA, Fairweather SJ, MacDermott-Opeskin HI, et al. (2021) The role of plasmalogens, Forssman lipids, and sphingolipid hydroxylation in modulating the biophysical properties of the epithelial plasma membrane. J Chem Phys 154: 095101. https://doi.org/10.1063/5.0040887
[41]
Marrink SJ, Risselada HJ, Yefimov S, et al. (2007) The MARTINI force field: coarse grained model for biomolecular simulations. J Phys Chem B 111: 7812-7824. https://doi.org/10.1021/jp071097f
[42]
de Jong DH, Singh G, Bennett WFD, et al. (2013) Improved parameters for the Martini coarse-grained protein force field. J Comp Theory Comput 9: 687-697. https://doi.org/10.1021/ct300646g
[43]
Abraham MJ, Murtola T, Schulz R, et al. (2015) GROMACS: high performance molecular simulations through multi-level parallelism from laptops to supercomputers. Software X 1: 19-25. https://doi.org/10.1016/j.softx.2015.06.001
[44]
Marrink SJ, de Vries AH, Mark AE (2004) Coarse grained model for semiquantitative lipid simulations. J Phys Chem B 108: 750-760. https://doi.org/10.1021/jp036508g
Parchekani J, Allahverdi A, Taghdir M, et al. (2022) Design and simulation of the liposomal model by using a coarse-grained molecular dynamics approach towards drug delivery goals. Sci Rep 12: 2371. https://doi.org/10.1038/s41598-022-06380-8
[50]
Pezeshkian W, König M, Wassenaar TA, et al. (2020) Backmapping triangulated surfaces to coarse-grained membrane models. Nat Commun 11: 2296. https://doi.org/10.1038/s41467-020-16094-y
Arnarez C, Uusitalo JJ, Masman MF, et al. (2015) Dry Martini, a coarse-grained force field for lipid membrane simulations with implicit solvent. J Chem Theory Comput 11: 260-275. https://doi.org/10.1021/ct500477k
[54]
Kučerka N, Nieh MP, Katsaras J (2011) Fluid phase lipid areas and bilayer thicknesses of commonly used phosphatidylcholines as a function of temperature. Biochim Biophys Acta 1808: 2761-2771. https://doi.org/10.1016/j.bbamem.2011.07.022
[55]
Petrache HI, Tristram-Nagle S, Gawrisch K, et al. (2004) Structure and fluctuations of charged phosphatidylserine bilayers in the absence of salt. Biophys J 86: 1574-1586. https://doi.org/10.1016/S0006-3495(04)74225-3
[56]
Petrache HI, Dood SW, Brown MF (2000) Area per lipid and acyl length distributions in fluid phosphatidylcholines determined by 2H NMR spectroscopy. Biophys J 79: 3172-3192. https://doi.org/10.1016/S0006-3495(00)76551-9
[57]
Brown AR, Sachs JN (2014) Determining structural and mechanical properties from molecular dynamics simulations of lipid vesicles. J Chem Theory Comput 10: 4160-4168. https://doi.org/10.1021/ct500460u
[58]
Gurtovenko AA, Vattulainen I (2008) Effect of NaCl and KCl on phosphatidylcholine and phosphatidylethanolamine lipid membranes: insight from atomic-scale simulations for understanding salt-induced effects in the plasma membrane. J Phys Chem B 112: 1953-1962. https://doi.org/10.1021/jp0750708
[59]
Lui A, Qi X (2012) Molecular dynamics simulations of DOPC lipid bilayers: the effect of Lennard-Jones parameters of hydrocarbon chains. Comput Mol Biosci 2: 78-82. https://doi.org/10.4236/cmb.2012.23007
Jurkiewicz P, Cwiklik L, Vojtišková A, et al. (2012) Structure, dynamics, and hydration of POPC/POPS bilayers suspended in NaCl, KCl, and CsCl solutions. Biochim Biophys Acta 1818: 609-616. https://doi.org/10.1016/j.bbamem.2011.11.033
[62]
Hinton DP, Johnson CS (1993) Diffusion ordered 2D NMR spectroscopy of phospholipid vesicles: determination of vesicle size distributions. J Phys Chem 97: 9064-9072. https://doi.org/10.1021/j100137a038
[63]
Muhsin-Sharafaldine MR, Saunderson SC, Dunn AC, et al. (2016) Procoagulant and immunogenic properties of melanoma exosomes, microvesicles and apoptotic vesicles. Oncotarget 7: 56279-56294. https://doi.org/10.18632/oncotarget.10783
[64]
Pereno V, Carugo D, Bau L, et al. (2017) Electroformation of giant unilamellar vesicles on stainless steel electrodes. ACS Omega 2: 994-1002. https://doi.org/10.1021/acsomega.6b00395
[65]
Lin CM, Wu DT, Tsao HK, et al. (2012) Membrane properties of swollen vesicles: growth, rupture, and fusion. Soft Matter 8: 6139. https://doi.org/10.1039/c2sm25518a
[66]
Olsson U, Nakamura K, Kunieda H, et al. (1996) Normal and reverse vesicles with nonionic surfactant: solvent diffusion and permeability. Langmuir 12: 3045-3054. https://doi.org/10.1021/la9600560
[67]
Przybylo M, Sýkora J, Humpolíčová J, et al. (2006) Lipid diffusion in giant unilamellar vesicles is more than 2 times faster than in supported phospholipid bilayers under identical conditions. Langmuir 22: 9096-9099. https://doi.org/10.1021/la061934p
[68]
Macháň R, Hof M (2010) Lipid diffusion in planar membranes investigated by fluorescence correlation spectroscopy. Biochim Biophys Acta 1798: 1377-1391. https://doi.org/10.1016/j.bbamem.2010.02.014
[69]
Róg T, Murzyn K, Pasenkiewicz-Gierula M (2003) Molecular dynamics simulations of charged and neutral lipid bilayers: treatment of electrostatic interactions. Acta Biochim Pol 50: 789-798.
[70]
Wang Y, Markwick PR, de Oliveira CA, et al. (2011) Enhanced lipid diffusion and mixing in accelerated molecular dynamics. J Chem Theory Comput 7: 3199-3207. https://doi.org/10.1021/ct200430c
[71]
Akhunzada MJ, D'Autilia F, Chandramouli B, et al. (2019) Interplay between lipid lateral diffusion, dye concentration and membrane permeability unveiled by a combined spectroscopic and computational study of a model lipid bilayer. Sci Rep 9: 1508. https://doi.org/10.1038/s41598-018-37814-x
[72]
Almeida PF, Vaz WLC, Thompson TE (1992) Lateral diffusion in the liquid phases of dimyristoylphosphatidylcholine/cholesterol lipid bilayers: a free volume analysis. Biochemistry 31: 6739-6747. https://doi.org/10.1021/bi00144a013
[73]
Kalli AC, Rog T, Vattulainan I, et al. (2017) The integrin receptor in biologically relevant bilayers: insights from molecular dynamics simulations. J Membrane Biol 250: 337-351. https://doi.org/10.1007/s00232-016-9908-z
[74]
Lira RB, Steinkühler J, Knorr RL, et al. (2016) Posing for a picture: vesicle immobilization in agarose gel. Sci Rep 6: 25254. https://doi.org/10.1038/srep25254
[75]
Wang T, Ingram C, Weisshaar JC (2010) Model lipid bilayer with facile diffusion of lipids and integral membrane proteins. Langmuir 26: 11157-11164. https://doi.org/10.1021/la101046r
Stier A, Sackmann E (1973) Spin labels as enzyme substrates. Heterogeneous lipid distribution in liver microsomal membranes. Biochim Biophys Acta 311: 400-408. https://doi.org/10.1016/0005-2736(73)90320-9
Semrau S, Schmidt T (2009) Membrane heterogeneity–from lipid domains to curvature effects. Soft Matter 5: 3174-3186. https://doi.org/10.1039/b901587f
[80]
Gu RX, Baoukina S, Tieleman DP (2019) Cholesterol flip-flop in heterogeneous membranes. J Chem Theory Comput 15: 2064-2070. https://doi.org/10.1021/acs.jctc.8b00933
[81]
Ermilova I, Lyubartsev AP (2019) Cholesterol in phospholipid bilayers: positions and orientations inside membranes with different unsaturation degrees. Soft Matter 15: 78-93. https://doi.org/10.1039/c8sm01937a
[82]
Buwaneka P, Ralko A, Liu SL, et al. (2021) Evaluation of the available cholesterol concentration in the inner leaflet of the plasma membrane of mammalian cells. J Lipid Res 62: 100084. https://doi.org/10.1016/j.jlr.2021.100084
[83]
Mondal M, Mesmin B, Mukherjee S, et al. (2009) Sterols are mainly in the cytoplasmic leaflet of the plasma membrane and the endocytic recycling compartment in CHO cells. Mol Biol Cell 20: 581-588. https://doi.org/10.1091/mbc.e08-07-0785
[84]
Courtney KC, Pezeshkian W, Raghupathy R, et al. (2018) C24 sphingolipids govern the transbilayer asymmetry of cholesterol and lateral organization of model and live-cell plasma membranes. Cell Rep 24: 1037-1049. https://doi.org/10.1016/j.celrep.2018.06.104
[85]
Oh Y, Sung BJ (2018) Facilitated and non-Gaussian diffusion of cholesterol in liquid ordered phase bilayers depends on the flip-flop and spatial arrangement of cholesterol. J Phys Chem Lett 9: 6529-6535. https://doi.org/10.1021/acs.jpclett.8b02982
[86]
Martínez JM, Leandro Martínez L (2003) Packing optimization for automated generation of complex system's initial configurations for molecular dynamics and docking. J Comput Chem 24: 819-825. http://doi.org/10.1002/jcc.10216
[87]
Martínez L, Andrade R, Birgin EG, et al. (2009) PACKMOL: a package for building initial configurations for molecular dynamics simulations. J Comput Chem 30: 2157-2164. https://doi.org/10.1002/jcc.21224
[88]
Jefferys E, Sands ZA, Shi J, et al. (2015) Alchembed: A computational method for incorporating multiple proteins into complex lipid geometries. J Chem Theory Comput 11: 2743-2754. https://doi.org/10.1021/ct501111d
[89]
Krause MR, Regen SL (2014) The structural role of cholesterol in cell membranes: from condensed bilayers to lipid rafts. Acc Chem Res 47: 3512-3521. https://doi.org/10.1021/ar500260t
[90]
Levental I, Levental KR, Heberle FA (2020) Lipid rafts: controversies resolved, mysteries remain. Trends Cell Biol 30: 341-353. https://doi.org/10.1016/j.tcb.2020.01.009
[91]
Reddy AS, Warshaviak DT, Chachisvilis M (2012) Effect of membrane tension on the physical properties of DOPC lipid bilayer membrane. Biochim Biophys Acta 1818: 2271-2281. https://doi.org/10.1016/j.bbamem.2012.05.006
Chia-Wen Wang, Meng-Han Lin, Wolfgang B. Fischer. Cholesterol affected dynamics of lipids in tailor-made vesicles by ArcVes software during multi micro second coarse grained molecular dynamics simulations[J]. AIMS Biophysics, 2023, 10(4): 482-502. doi: 10.3934/biophy.2023027
Chia-Wen Wang, Meng-Han Lin, Wolfgang B. Fischer. Cholesterol affected dynamics of lipids in tailor-made vesicles by ArcVes software during multi micro second coarse grained molecular dynamics simulations[J]. AIMS Biophysics, 2023, 10(4): 482-502. doi: 10.3934/biophy.2023027
Figure 1. Arc-corrected Vesicle construction (ArcVes): The two entry windows to generate tailor-made vesicle by placing the lipid molecule and proteins as well as ‘customizing’ the distributing of the lipid molecules in each leaflet
Figure 2. Coarse grained (CG) models of mixed vesicles with a setup-radius of 10 nm. (A) Upper panel: vesicles of POPC/cholesterol (POPC/CHOL) generated by ArcVes (left) and at the end of the simulation (right, 10 µs) with POPC represented in grey (outer leaflet) and light blue spheres (inner leaflet). Lower panel: vesicles of DOPC/cholesterol (DOPC/CHOL) with DOPC represented in green and yellow for outer and inner leaflet, respectively. Cholesterol is represented by blue (outer leaflet) and red (inner leaflet) spheres with the ROH spheres highlighted in magenta for outer and black for inner leaflet to represent the head group of cholesterol. The vesicles present insight into the interior by omitting half of the molecules of the outer leaflet and one quarter of the molecules of the inner leaflets. (B) The same mixed vesicles showing the lipids separately for respective POPC and DOPC (‘PO4’-spheres in blue for outer and red for inner leaflet) and the respective cholesterol molecules (‘ROH’-spheres in blue for outer and red for inner leaflet)
Figure 3. Time dependent calculate properties from the simulations. (A) Radius [nm], (B) thickness [nm], and (C) number of lipids [count × 102]. The left graphs represent the mixed vesicle POPC/CHOL, and the right graphs present the DOPC/CHOL. Thickness values for POPC and DOPC are shown in green line and cholesterol in red line. Standard deviations (SD [nm]) are shown in grey area representation in the radius plot and for POPC and DOPC in thickness plot. The SD for cholesterol is shown by a black line in the thickness plots. The value for the ArcVes setup-radius is marked as a light grey dashed line. In the plots for the number of lipids, values for POPC, DOPC and cholesterol are all shown in blue and red for outer (o) and inner leaflets (i), respectively
Figure 4. 2D projection of the density [count/nm2] of the lipids in the mixed vesicles POPC/CHOL (left) and DOPC/CHOL (right) calculated from a 10 µs MD simulation. (A) Density of both leaflets. (B) Density of the outer (upper row, o) ad inner leaflet (lower row, i). The calculation is done by flattening the entire vesicles resulting in axis ranging from 0 to 180 ° and to 360 °. The color coding goes from low to high density as black < dark blue < light blue < green < yellow < red < white