Stochastic Variational Methods in Generalized Hidden Semi-Markov Models to Characterize Functionality in Random Heteropolymers

Protein-mimicry with synthetic polymers has been actively pursued for decades to meet material demands within the environmental, energy and life sciences. Recently, Panganiban et al. (2018) and Jiang et al. (2020) reported a promising protein-mimic polymer system on the basis of random heterogeneous polymers (RHP) comprising monomers with distinct chemical properties. Through wet lab experiments, they found long-term screening for optimizing the monomer choices and monomer ratios can lead to some RHPs with protein-like functions. The RHP systems provide an attractive platform for developing polymer-based biofunctional materials.

In Jiang et al. (2020), several datasets of four-monomer based RHPs were synthesized, which behave like natural protein-based proton channels. The four methacrylate-based monomers were selected by chemical diversity, i.e., methyl methacrylate (MMA), 2-ethylhexyl methacrylate (EHMA), oligo (ethylene glycol) methacrylate (OEGMA), and 3-sulfopropyl methacrylate potassium salt (SPMA). MMA and EHMA are the hydrophobic monomers, while OEGMA and SPMA are the hydrophilic monomers. RAFT polymerization was used to control the ratio of four monomers in an RHP dataset, in order to tailor the overall hydrophobicity. Different RHP datasets with varying monomer ratios were generated, and some of them were found capable of being inserted into lipid bilayers and promoting transmembrane proton transport. However, there is a difference from natural proteins that make up the channels. A functional natural protein has certain specific amino acid sequence(s). But the functional performance of each RHP dataset is determined as a whole, and chains in an RHP dataset, which are synthesized under a fixed monomer ratio, are largely different from each other. The sequence heterogeneity in RHP brings challenges to directly linking individual monomer sequence to the overall performance of a polymer dataset. Study of the sequence heterogeneity under various monomer composition ratios would shed light on the understanding of RHPs’ functional adaptability, enabling their uniform behavior in various environments. More importantly, it offers promising potential for further RHP design in a predictable manner.

As RHPs’ local hydrophobicity levels relate to short-range interactions, each RHP chain can be studied as a concatenation of segments of various cumulative hydrophobicity. It is believed that varied segment lengths and ending positions substantially affect RHPs’ chemical properties. Therefore, segmentation can reveal valuable latent structure information on RHP heterogeneity and functionalities when individual chain experimentation is unavailable. The differences in cumulative hydrophobicity among segments introduce the segment heterogeneity. To analyze such heterogeneity, the segments are assigned to three states: the hydrophilic state, $S1$, for the segments rich in hydrophilic monomer OEGMA and the negatively charged SPMA; the hydrophobic state, $S3$, for those comprising mostly hydrophobic monomers MMA and EHMA; the intermediate state, $S2$, for ones that cannot be simply assigned as hydrophobic or hydrophilic (amphiphilic segments). The $S3$ segments have the ability to span across the bilayer region and provide the proton transport pathway. The $S1$ segments prefer the water environment. The $S2$ segments are believed to have short residence time in both bilayer and water regions.

Figure 1 displays the schematic illustration of RHP chains in a lipid bilayer. Note however, information on segmentation and spatial structure is not available in practice, and we only have access to the linear sequence of four monomers in each RHP dataset.

Determining how RHPs’ performance relates to segment heterogeneity and sequence heterogeneity, rather than to a specific realization of polymer sequence, provides a promising strategy for developing biofunctional materials. The two types of heterogeneity correspond to two key factors in RHP design: (1) the segmentation of local hydrophobicity that regulates short-range interactions, and (2) the choice of overall monomer composition ratio that modulates the solubility. So far, domain knowledge is extremely limited in the field. There is also a lack of statistical study that could bridge the experimental restrictions and RHPs’ predictable performance. To exploit the statistical design of RHP and to aid its synthesis, we utilized the probabilistic graphical model framework and developed a generalized hidden semi-Markov model (GHSMM-RHP). Based on this model, we predicted the segments of varying hydrophobicity within RHP chains to address the first key factor above. To address the second factor, we statistically quantified the heterogeneity of chain distributions among different RHP datasets. GHSMM-RHP helps to extract the hidden phenomena that are otherwise hard to quantify, and utilizes efficient approximate inference in place of expensive exact computation. That is, we aim to address a joint estimation-prediction problem in RHP systems under the GHSMM-RHP model.

Specifically, GHSMM-RHP adopts the residential time design, which has fewer parameters and lower time complexity than generic hidden semi-Markov model (Yu, 2010). GHSMM-RHP incorporates additional layer of nodes to allow handling segment constraints such as the monomer positions, as outlined in Section 3. Based on the model, we developed stochastic variational methods in closed-form under both frequentist and Bayesian settings. To keep the monomer dependencies within an RHP chain, we derived message-passing formulas on the graphical model of conjugate exponential families (Maathuis et al., 2018). GHSMM-RHP differs from existing efforts in its novel design that incorporates specific constraints for large RHP systems, and provides a scalable graphical model with efficient algorithms that meets the inference demand in RHP applications.

In addition, we studied the GHSMM-RHP model through a simulation analysis, and compared estimates from the frequentist and Bayesian stochastic variational methods. In particular, we implemented stochastic variational expectation-maximization (SVEM) and stochastic variational Bayesian (SVB), and empirically addressed their estimation/prediction accuracy and computational (in)stability against hyper-parameters. We found that both methods achieved low statistical error, while SVEM exhibited greater stability against random initializations. Our work makes the first attempt to bring the two stochastic variational methods, which were studied separately before, to the same nonconvex regime. It also contributes to the current literature on the computational error of stochastic variational methods.

We further evaluated the GHSMM-RHP model using real datasets. Jiang et al. (2020) synthesized seven RHP datasets with varied compositional ratios and tested their performance. Among them, RHP with a compositional ratio of 50(MMA) : 20(EHMA) : 25(OEGMA) : 5(SPMA) was the most functionally active, which we refer to as RHP50/20. We used RHP50/20 as a reference and compared other RHPs to RHP50/20 based on the GHSMM-RHP model. The RHP datasets identified by our model as having high similarity to RHP50/20 all produced high proton transport efficiency in the lab experiment, which suggests that the GHSMM-RHP model has the potential to aid in predictable RHP performance. By bridging the gap between costly laboratory assessments and less expensive computational models, we hope to provide guidance for designing improved RHP syntheses.

The rest of the paper is organized as follows. In Section 2 we review the background of methods related to the GHSMM-RHP model. In Section 3 we describe the model that incorporates the constraints of RHP systems. In Section 4 different stochastic variational methods are derived. We implemented the algorithms and compared the results in Section 5. Section 6 concludes the work and discusses avenues for further research.

As described previously, each RHP dataset contains thousands of sequences randomly generated under a fixed monomer ratio. Note that we only have access to the one-dimensional primary structure, that is, the linear sequence of four monomers. We start by discussing a few sequence constraints, which are the same as those stated in Jiang et al. (2020). These constraints describe the underlying biology of functional RHPs: (1) We define $S3$ segments to be composed of hydrophobic monomers (i.e., MMA, EHMA) and no more than one OEGMA. If an OEGMA appears in an S3 segment, it must be at least two monomers away from both ends of the $S3$ segment. (2) $S2$ and $S3$ are separated by hydrophilic segments $S1$. (3) The negatively charged hydrophilic monomer SPMA appears only in $S1$.

The hydrophobic monomers in $S3$ segments enable the insertion of the chain into lipid bilayers, and the appearance of a single OEGMA doesn’t compromise the overall activity if it is in the middle of an $S3$ segment. The length of the $S3$ segments depends on various factors such as chain conformation and membrane flexibility, but it should be sufficiently long to span the core of the lipid bilayer. The $S1$ segments have a high preference for water (high polarity of the monomer), and appear in between $S2$ and $S3$.

We adapted the right-censored hidden semi-Markov model in (Yu, 2010) to model the RHP chain while incorporating the above constraints. Here, the hidden state duration is a random variable following a categorical distribution over a predefined range. The unobserved states correspond to three types of segments: hydrophilic state $S1$, intermediate state $S2$, and hydrophobic state $S3$. We denote a random RHP chain as $(y_{1},...,y_{K})$. At each position $k\in 1,...,K$, $y_{k}$ takes values from the four monomers, i.e., MMA, OEGMA, EHMA, SPMA. $(z_{k},\tau_{k},\iota_{k})$ are the latent variables. $z_{k}$ is the hidden segment state, $\tau_{k}$ is the remaining duration of the current hidden state, and $\iota_{k}$ is the location of an emitted OEGMA. Their distributions are described below.

$\bullet\hskip 2.84526ptz_{1}$: The model starts with the initial state variable $z_{1}$, taking value in the state space $\mathcal{S}_{z}=\{S1,S2,S3\}$. $p(z_{1}|\pi)$ follows a categorical distribution with parameters $\pi=(\pi_{S1},\pi_{S2},\pi_{S3})$:

$\bullet\hskip 2.84526pt(\tau_{1},\iota_{1})$: Latent variables $(\tau_{1},\iota_{1})$ depend on $z_{1}$. (1) If $z_{1}=S3$, then the duration $\tau_{1}$ takes values between $\underline{D}$ and $\overline{D}$, the predefined shortest and longest segment length. The OEGMA location $\iota_{1}$ takes values between 3 and $\tau_{1}-2$, meaning the OEGMA shall be at least two monomers away from both ends of the $S3$ segment. $\iota_{1}$ may also be 0, meaning there is no OEGMA in the current $S3$ segment. This support is denoted as $\mathcal{S}^{S3}_{\tau,\iota}=\{(d,l)|(d,l)\in\{\underline{D},...,\overline{D}\}\times\{0,3,...,\overline{D}-2\},l\leq d-2\}$. The conditional probability $p(\tau_{1},\iota_{1}|z_{1},\rho^{S3})$ follows a categorical distribution with parameters $\rho^{S3}=\{\rho^{S3}_{d,l}|(d,l)\in\mathcal{S}^{S3}_{\tau,\iota}\}$. (2) If $z_{1}=S1$ or $z_{1}=S2$, then the duration $\tau_{1}$ takes value between $1$ and $\overline{D}$. There is no constraint on OEGMA location for $S1$ and $S2$ segments, so $\iota_{1}$ can take any value without affecting the likelihood function. We set it to 0 for convenience. This support is denoted as $\mathcal{S}^{S1}_{\tau,\iota}=\mathcal{S}^{S2}_{\tau,\iota}=\{(d,l)|(d,l)\in\{1,...,\overline{D}\}\times\{0\}\}$. The conditional probability $p(\tau_{1},\iota_{1}|z_{1},\rho^{S1})$ follows a categorical distribution with parameters $\rho^{S1}=\{\rho^{S1}_{d,l}|(d,l)\in\mathcal{S}^{S1}_{\tau,\iota}\}$. $p(\tau_{1},\iota_{1}|z_{1},\rho^{S2})$ and $\rho^{S2}$ is defined similarly.

In summary, the conditional probability $p(\tau_{1},\iota_{1}|z_{1},\rho)$ follows a categorical distribution with parameters $\rho=(\rho^{S1},\rho^{S2},\rho^{S3})$. Three sub-vectors $\rho^{S1},\rho^{S2},\rho^{S3}$ correspond to three supports $\mathcal{S}^{S1}_{\tau,\iota},\mathcal{S}^{S2}_{\tau,\iota},\mathcal{S}^{S3}_{\tau,\iota}$ respectively. This gives:

$\bullet\hskip 2.84526pt(z_{k},\tau_{k},\iota_{k})$: For positions $k=2,...,K$, latent variables $(z_{k},\tau_{k},\iota_{k})$ evolve as a time-homogeneous Markov chain, taking values in the same joint space as the one described for $(z_{1},\tau_{1},\iota_{1})$. The graphical model diagram is illustrated in Figure 2. The transition distribution with parameters $\alpha,\rho$ is defined as $\zeta_{(z_{k-1},\tau_{k-1},\iota_{k-1})\to(z_{k},\tau_{k},\iota_{k})}$. When $(z_{k},\tau_{k},\iota_{k})=(j,d^{\prime},l^{\prime})$ and $(z_{k-1},\tau_{k-1},\iota_{k-1})=(i,d,l)$, the transition kernel $\zeta$ is:

As the position $k$ proceeds we note the following. (1) If $d>1$, there is no transition on hidden segment state, i.e., $z_{k}=j=i$. The remaining duration $\tau_{k-1}=d$ of hidden segment state $z_{k-1}=i$ will deterministically decrease by 1 until it reaches 1, i.e., $\tau_{k}=d-1$, $\tau_{k+1}=d-2$, et cetera. (2) If $d=1$, the RHP will transit to a different hidden segment state $z_{k}=j\in\mathcal{S}_{z}\setminus\{i\}$ with probability $\alpha_{i,j}$, followed by a new pair $(\tau_{k},\iota_{k})=(d^{\prime},l^{\prime})$ generated according to the categorical probability parameter $\rho^{j}$. We disallow transitions between S2 and S3 due to biological constraints, i.e., $\alpha_{S2,S3}=\alpha_{S3,S2}=0$. Thus, we have:

$\bullet\hskip 2.84526pty_{k}$: The observed random variable $y_{k}$ is conditionally independent of all the other variables given $(z_{k},\tau_{k},\iota_{k})$. $y_{k}$ takes values in the space of four monomers $\mathcal{S}_{y}=\{MMA,OEGMA,EHMA,SPMA\}$. The emission distribution $p(y_{k}|z_{k},\tau_{k},\iota_{k},\beta)$ follows a categorical distribution with parameters $\beta=(\beta_{S1},\beta_{S2},\beta_{S3},\beta_{S3^{\prime}})$. The four sub-vectors of $\beta$ are used conditional on different values of $(z_{k},\tau_{k},\iota_{k})$, as shown in Table 1.

Specifically, $\beta_{S1}$ is for categorical distribution over $\mathcal{S}_{y}$ when $z_{k}=S1$. $\beta_{S2}$ is for categorical distribution over $\mathcal{S}_{y}\setminus\{SPMA\}$ when $z_{k}=S2$, because we assume the SPMA monomer appears only in $S1$ segments. $\beta_{S3}$ is for categorical distribution over $\mathcal{S}_{y}\setminus\{OEGMA,SPMA\}$ if $z_{k}=S3,\tau_{k}\neq\iota_{k}$ or $z_{k}=S3,\iota_{k}=0$. It is used when the current position $k$ is within an $S3$ segment but no OEGMA is emitted, either because $\tau_{k}$ hasn’t reached the designated OEGMA location $\iota_{k}$ or because the current $S3$ segment is made up of purely MMA and EHMA. $\beta_{S3^{\prime}}$ is for categorical distribution over $\{OEGMA\}$ if $z_{k}=S3,\tau_{k}=\iota_{k}$. It degenerates to a deterministic value. It is used when the duration variable $\tau_{k}$ decreases to the designated OEGMA location $\iota_{k}$, indicating that the only OEGMA should be emitted in the current $S3$ segment. Thus we have:

This probabilistic model incorporates the aforementioned constraints on RHP systems, and allows for explicit modeling of the hidden states’ duration. The complete likelihood function for one RHP chain is:

In the last line of Equation LABEL:eq:complete_likelihood, we plug in Equations 2 to LABEL:eq:emission and simplify Equation LABEL:eq:complete_likelihood further to an exponential family. Here, ${\theta}=(\pi,\alpha,\beta,\rho)$. ${\eta}({\theta})=\ln\theta$ is the canonical parameter vector and $\ln$ is applied element-wise. We denote $\bm{y}=y_{1:K}=(y_{1},...,y_{K})$, with $\bm{z,\tau,\iota}$ similarly defined. ${\phi}(\bm{z,\tau,\iota,y})$ is the vector of complete-data sufficient statistics. All sequences are i.i.d. given the parameter ${\theta}$, but for clarity we have suppressed indices for i.i.d. replicas of RHP sequences so far. Denote $\bm{y}^{(s)}=y_{1:K}^{(s)}=(y_{1}^{(s)},...,y_{K}^{(s)})$ as observations of the $s$-th sequence, $s=1,...,N$. Denote $\bm{y}^{(\mathcal{I})}=\{\bm{y}^{(s)}\}_{s\in\mathcal{I}}$ as all RHP sequences in the index set $\mathcal{I}\subseteq\{1,...,N\}$, i.e., $\bm{y}^{(1:N)}=\{\bm{y}^{(s)}\}_{s=1}^{N}$. $\bm{z}^{(s)},\bm{\tau}^{(s)},\bm{\iota}^{(s)}$ and $\phi^{(s)}=\phi(\bm{z}^{(s)},\bm{\tau}^{(s)},\bm{\iota}^{(s)},\bm{y}^{(s)})$ are defined similarly.

Variational methods have been used for approximate inference under different scenarios when exact solutions are intractable or too expensive to achieve (Bishop, 2006; Wainwright and Jordan, 2008). To develop more efficient algorithms in modern large-scale learning problems, stochastic optimizations are introduced to bypass the need to iterate through the entire datasets in each iteration (Robbins and Monro, 1951). Combining variational methods and stochastic optimizations helps alleviate the computational complexity of estimations in the GHSMM-RHP model.

In particular, in the Bayesian setting where parameters $\theta$ are viewed as random variables, we coupled stochastic natural gradient descent (Amari, 1998) with variational Bayes following the work of Hoffman et al. (2013). In the frequentist setting where parameters $\theta$ are viewed as fixed but unknown, we coupled the stochastic Frank-Wolfe algorithm (Reddi et al., 2016) with variational expectation-maximization.

We implemented SVB and SVEM, and compared their performance using a comprehensive simulation study and a real data study from Jiang et al. (2020).