chain_model.cpp
1 #include <cmath>
2 #include <stdlib.h>
3 #include "chain_model.h"
4 #include "parsing.h"
5 
6 namespace BioLCCC
7 {
8 
9 std::vector<double> calculateBoltzmannFactorProfile(
10  const std::vector<double> &effectiveEnergyProfile)
11 {
12  std::vector<double> boltzmannFactorProfile;
13 
14  for (std::vector<double>::const_iterator energy =
15  effectiveEnergyProfile.begin();
16  energy != effectiveEnergyProfile.end();
17  ++energy)
18  {
19  boltzmannFactorProfile.push_back(exp(*energy));
20  }
21 
22  return boltzmannFactorProfile;
23 }
24 
26  const std::vector<ChemicalGroup> &parsedSequence,
27  const double secondSolventConcentration,
28  const ChemicalBasis &chemBasis,
29  const double columnPoreSize,
30  const double columnRelativeStrength,
31  const double temperature
32  ) throw (BioLCCCException)
33 {
34  // At first, we need to convert the energy profile to a profile of
35  // distribution probabilities. Probability = exp(E_effective),
36  // where E_effective = E_of_residue - Eab,
37  // and Eab is an energy of binding for a binary solvent
38  // and Eab = Ea + ln ( 1 + Nb + Nb*exp (Ea - Eb) )
39  // also corrections of energies due to temperature (energies in exponents
40  // are scaled to the new temperature) and the relative strength
41  // are introduced.
42  std::vector<std::vector<double> > boltzmannFactorProfiles;
43  for (std::vector<double>::const_iterator adsorptionLayerFactor =
44  chemBasis.adsorptionLayerFactors().begin();
45  adsorptionLayerFactor != chemBasis.adsorptionLayerFactors().end();
46  ++adsorptionLayerFactor)
47  {
48  boltzmannFactorProfiles.push_back(
52  parsedSequence,
53  chemBasis,
54  secondSolventConcentration,
55  columnRelativeStrength * (*adsorptionLayerFactor),
56  temperature),
57  chemBasis.monomerLength(),
58  chemBasis.kuhnLength())));
59  }
60 
61  // The size of the lattice must be greater than
62  // (number of adsorbing layers) * 2.
63  // double round (double x) {return floor(x+0.5);}
64  unsigned int latticeSize =
65  floor(columnPoreSize / chemBasis.kuhnLength() + 0.5);
66 
67  // If we want to neglect the partially desorbed states, we need to insert
68  // two impenetrable layers right after the near-wall ones.
69  if (chemBasis.neglectPartiallyDesorbedStates())
70  {
71  boltzmannFactorProfiles.push_back(
72  std::vector<double>(
73  boltzmannFactorProfiles.back().size(), 0.0));
74  latticeSize += 2;
75  }
76 
77  if (latticeSize < boltzmannFactorProfiles.size() * 2)
78  {
79  throw BioLCCCException(
80  "The pore size is too small for the given number of adsorbing "
81  "layers.");
82  }
83 
84  // The density vector correspond to a probability of n-th residue to be in
85  // a certain layer between pore walls.
86  double *density;
87 
88  // The transition matrix used to calculate a density vector of n-th
89  // particle from a density vector of (n-1)-th particle.
90  double *transitionMatrix;
91 
92  // The density buffer vector is used during matrix calculations.
93  double *densityBuffer;
94 
95  // Memory managment.
96  try
97  {
98  density = new double[latticeSize];
99  densityBuffer = new double[latticeSize];
100  transitionMatrix = new double[latticeSize*latticeSize];
101  }
102  catch (...)
103  {
104  throw BioLCCCException("Cannot allocate memory for calculations");
105  }
106 
107  // Constructing a density vector for the first amino acid residue.
108  // A density is distributed uniformly over all non-adsorbing layers of
109  // the lattice.
110  // The density in adsorbing layers is multiplied by a Boltzmann factor of
111  // the first segment.
112 
113  for (unsigned int i = 0; i < latticeSize; i++)
114  {
115  density[i] = 1.0;
116  }
117 
118  for (unsigned int i = 0; i < boltzmannFactorProfiles.size(); ++i)
119  {
120  density[i] = boltzmannFactorProfiles[i][0];
121  density[latticeSize - i - 1] = boltzmannFactorProfiles[i][0];
122  }
123 
124  // Debugging facilities.
125  //for (unsigned int i = 0; i < latticeSize; i++)
126  //{
127  // std::cout << density[i] << " ";
128  //}
129  //std::cout << std::endl;
130  //std::cout << std::endl;
131 
132  // Than we construct a basis for the transition matrix. The basis is
133  // a diagonal matrix with 4.0/6.0 on the main diagonal and 1.0/6.0 on
134  // the side diagonals.
135 
136  // Filling the matrix.
137  for (unsigned int i = 0; i < latticeSize; i++)
138  {
139  for (unsigned int j = 0; j < latticeSize; j++)
140  {
141  switch ((int)abs((int)(j - i)))
142  {
143  case 0:
144  {
145  transitionMatrix[j + latticeSize * i] = 4.0/6.0;
146  break;
147  }
148  case 1:
149  {
150  transitionMatrix[j + latticeSize * i] = 1.0/6.0;
151  break;
152  }
153  default:
154  transitionMatrix[j + latticeSize * i] = 0.0;
155  }
156  }
157  }
158 
159  // On each step we obtain the density vector for the n-th segment
160  // multiplying the transition matrix and the density vector of the
161  // (n-1)th residue.
162  // The multiplication starts from the second segment.
163  for (unsigned int segmentIndex = 1;
164  segmentIndex < boltzmannFactorProfiles[0].size();
165  ++segmentIndex)
166  {
167  // Filling the matrix elements that describe the adsorption.
168  for (unsigned int layerIndex = 0;
169  layerIndex < boltzmannFactorProfiles.size();
170  ++layerIndex)
171  {
172  int indexShift = layerIndex * ( latticeSize + 1 );
173  double boltzmannFactor =
174  boltzmannFactorProfiles[layerIndex][segmentIndex];
175 
176  transitionMatrix[indexShift + 0] = 4.0/6.0 * boltzmannFactor;
177  transitionMatrix[indexShift + 1] = 1.0/6.0 * boltzmannFactor;
178  transitionMatrix[latticeSize*latticeSize - 1 - indexShift - 0] =
179  4.0 / 6.0 * boltzmannFactor;
180  transitionMatrix[latticeSize*latticeSize - 1 - indexShift - 1] =
181  1.0 / 6.0 * boltzmannFactor;
182 
183  // A segment may enter the second and further adsorption layers from
184  // the inner layer (i.e. the layer lying closer to the walls).
185  if (layerIndex > 0)
186  {
187  transitionMatrix[indexShift - 1] = 1.0/6.0 * boltzmannFactor;
188  transitionMatrix[latticeSize*latticeSize - 1 - indexShift + 1] =
189  1.0 / 6.0 * boltzmannFactor;
190  }
191  }
192 
193  // Zeroing the calculation buffer.
194  for (unsigned int i = 0; i < latticeSize; i++)
195  {
196  densityBuffer[i] = 0.0;
197  }
198 
199  // Multiplying the transition matrix by the density vector. The result
200  // is stored in the buffer vector.
201  for (unsigned int i = 0; i < latticeSize; i++)
202  {
203  for (unsigned int j = 0; j < latticeSize; j++)
204  {
205  densityBuffer[i] = densityBuffer[i] + density[j] *
206  transitionMatrix[j + i * latticeSize];
207  }
208  }
209 
210  // Transferring the results from the density vector.
211  for (unsigned int i = 0; i < latticeSize; i++)
212  {
213  density[i] = densityBuffer[i];
214  }
215  }
216 
217  // Kd is a sum of elements of the density vector, normalized to the size
218  // of the lattice.
219  double Kd=0;
220  for (unsigned int i=0; i < latticeSize; i++)
221  {
222  Kd += density[i];
223  }
224 
225  if (chemBasis.neglectPartiallyDesorbedStates())
226  {
227  Kd = Kd / (double)(latticeSize - 2);
228  }
229  else
230  {
231  Kd = Kd / (double)(latticeSize);
232  }
233 
234  // Cleaning memory.
235  try
236  {
237  delete[] density;
238  delete[] densityBuffer;
239  delete[] transitionMatrix;
240  }
241  catch (...)
242  {
243  throw BioLCCCException("Cannot allocate memory for calculations");
244  }
245 
246  return Kd;
247 }
248 
249 }
An instance of ChemicalBasis contains a set of BioLCCC constants.
Definition: chemicalbasis.h:86
Apart from classes, BioLCCC contains calculation methods and constants.
Definition: auxiliary.h:4
std::vector< double > calculateMonomerEnergyProfile(const std::vector< ChemicalGroup > &parsedSequence, const ChemicalBasis &chemBasis, const double secondSolventConcentration, const double columnRelativeStrength, const double temperature)
Calculates the effective energy profile of monomers of the polymer chain.
Definition: parsing.cpp:11
Base class for all BioLCCC exceptions. Can be used by itself.
double calculateKdChain(const std::vector< ChemicalGroup > &parsedSequence, const double secondSolventConcentration, const ChemicalBasis &chemBasis, const double columnPoreSize, const double columnRelativeStrength, const double temperature)
Calculates coefficient of distribution of a polymer using the chain model.
Definition: chain_model.cpp:25
std::vector< double > calculateSegmentEnergyProfile(const std::vector< double > &monomerEnergyProfile, const double monomerLength, const double kuhnLength)
Calculates the effective energy profile of segments of the polymer chain.
Definition: parsing.cpp:87
std::vector< double > calculateBoltzmannFactorProfile(const std::vector< double > &effectiveEnergyProfile)
Converts energy profile of a peptide/protein to a profile of probabilities.
Definition: chain_model.cpp:9