biolccc.cpp
1 #include <vector>
2 #include <cmath>
3 #include <algorithm>
4 #include "biolccc.h"
5 
6 namespace BioLCCC
7 {
8 
9 // Auxiliary functions that shouldn't be exposed to a user.
10 namespace
11 {
12 double calculateKd(const std::vector<ChemicalGroup> &parsedSequence,
13  const double secondSolventConcentration,
14  const ChemicalBasis &chemBasis,
15  const double columnPoreSize,
16  const double columnRelativeStrength,
17  const double temperature) throw(BioLCCCException)
18 {
19  // assymetricCalculations shows whether the Kd for the reversed molecule
20  // will differ. It happens when a molecule cannot be divided into an integer
21  // number of Kuhn segments.
22  bool assymetricCalculations =
23  (fmod(chemBasis.monomerLength() * parsedSequence.size(),
24  chemBasis.kuhnLength()) != 0);
25  // Choosing the appropriate polymerModel.
26  if (chemBasis.polymerModel()==CHAIN)
27  {
28  double Kd = calculateKdChain(parsedSequence,
29  secondSolventConcentration,
30  chemBasis, columnPoreSize,
31  columnRelativeStrength, temperature);
32  if (assymetricCalculations)
33  {
34  std::vector<ChemicalGroup> revParsedSequence =
35  parsedSequence;
36  std::reverse(revParsedSequence.begin(),
37  revParsedSequence.end());
38  Kd = (Kd + calculateKdChain(revParsedSequence,
39  secondSolventConcentration,
40  chemBasis,
41  columnPoreSize,
42  columnRelativeStrength,
43  temperature)) / 2.0 ;
44  }
45  return Kd;
46  }
47  else if (chemBasis.polymerModel() == ROD)
48  {
49  double Kd = calculateKdRod(parsedSequence,
50  secondSolventConcentration,
51  chemBasis, columnPoreSize,
52  columnRelativeStrength, temperature);
53  if (assymetricCalculations)
54  {
55  std::vector<ChemicalGroup> revParsedSequence =
56  parsedSequence;
57  std::reverse(revParsedSequence.begin(),
58  revParsedSequence.end());
59  Kd = (Kd + calculateKdRod(revParsedSequence,
60  secondSolventConcentration,
61  chemBasis,
62  columnPoreSize,
63  columnRelativeStrength,
64  temperature)) / 2.0 ;
65  }
66  return Kd;
67  }
68  else
69  {
70  throw BioLCCCException("Model error.");
71  }
72 }
73 
74 class KdCalculator
75 {
76 public:
77  KdCalculator(const std::vector<ChemicalGroup> &parsedSequence,
78  const ChemicalBasis &chemBasis,
79  const double columnPoreSize,
80  const double columnRelativeStrength,
81  const double temperature,
82  const int numInterpolationPoints) :
83  mParsedSequence(parsedSequence),
84  mChemicalBasis(chemBasis),
85  mColumnPoreSize(columnPoreSize),
86  mColumnRelativeStrength(columnRelativeStrength),
87  mTemperature(temperature),
88  mN(numInterpolationPoints)
89  {
90  if (mN > 0)
91  {
92  mSecondSolventConcentrations = new double[mN];
93  mLogKds = new double[mN];
94  // The number of extra points in the terminal segments.
95  // This points significantly increase the accuracy of spline
96  // interpolation.
97  int NETP = 1;
98  for (int i=0; i < mN; i++)
99  {
100  if (i <= NETP)
101  {
102  mSecondSolventConcentrations[i] =
103  i * 100.0 / (mN - 2.0 * NETP - 1.0) / (NETP + 1.0);
104  }
105  else if (i > (mN - NETP - 2))
106  {
107  mSecondSolventConcentrations[i] =
108  ((mN - 2.0 * NETP - 2.0)
109  + (i - mN + NETP + 2.0) / (NETP + 1.0))
110  * 100.0 / (mN - 2.0 * NETP - 1.0);
111  }
112  else
113  {
114  mSecondSolventConcentrations[i] =
115  (i - NETP) * 100.0 / (mN - 2.0 * NETP - 1.0);
116  }
117  mLogKds[i] = log(calculateKd(mParsedSequence,
118  mSecondSolventConcentrations[i],
119  mChemicalBasis,
120  mColumnPoreSize,
121  mColumnRelativeStrength,
122  mTemperature));
123  }
124 
125  mSecondDers = new double[mN];
126  fitSpline(mSecondSolventConcentrations, mLogKds,
127  mN, mSecondDers);
128  }
129  }
130 
131  ~KdCalculator()
132  {
133  if (mN > 0)
134  {
135  delete[] mSecondSolventConcentrations;
136  delete[] mLogKds;
137  delete[] mSecondDers;
138  }
139  }
140 
141  double operator()(double secondSolventConcentration)
142  throw (BioLCCCException)
143  {
144  if (mN == 0)
145  {
146  return calculateKd(mParsedSequence,
147  secondSolventConcentration,
148  mChemicalBasis,
149  mColumnPoreSize,
150  mColumnRelativeStrength,
151  mTemperature);
152  }
153  else
154  {
155  return exp(calculateSpline(mSecondSolventConcentrations, mLogKds,
156  mSecondDers, mN,
157  secondSolventConcentration));
158  }
159  }
160 
161 private:
162  const std::vector<ChemicalGroup> &mParsedSequence;
163  const ChemicalBasis & mChemicalBasis;
164  const double mColumnPoreSize;
165  const double mColumnRelativeStrength;
166  const double mTemperature;
167  const int mN;
168  double * mSecondSolventConcentrations;
169  double * mLogKds;
170  double * mSecondDers;
171 };
172 
173 double calculateRT(const std::vector<ChemicalGroup> &parsedSequence,
174  const ChemicalBasis &chemBasis,
175  const ChromoConditions &conditions,
176  const int numInterpolationPoints,
177  const bool continueGradient,
178  const bool backwardCompatibility
179  ) throw(BioLCCCException)
180 {
181  if (numInterpolationPoints < 0)
182  {
183  throw BioLCCCException(
184  "The number of interpolation points must be non-negative.");
185  }
186 
187  KdCalculator kdCalculator(parsedSequence, chemBasis,
188  conditions.columnPoreSize(),
189  conditions.columnRelativeStrength(),
190  conditions.temperature(),
191  numInterpolationPoints);
192 
193  double RT = 0.0;
194  // Use simplified expression for isocratic elution.
195  if (conditions.SSConcentrations().size() == 1 || kdCalculator(conditions.SSConcentrations()[0]) <= 1.)
196  {
197  RT = kdCalculator(conditions.SSConcentrations()[0])
198  * conditions.columnPoreVolume() / conditions.flowRate();
199  }
200  else
201  {
202  // The part of a column passed by molecules. When it exceeds 1.0,
203  // the analyte elutes from the column.
204  double S = 0.0;
205  double dS = 0.0;
206  int j = 0;
207  double currentSSConcentration = 0.0;
208  while (S < 1.0)
209  {
210  j++;
211  if (j < conditions.SSConcentrations().size())
212  {
213  currentSSConcentration = conditions.SSConcentrations()[j];
214  }
215  else
216  {
217  // If continue gradient then use the slope of the last section.
218  if (continueGradient)
219  {
220  currentSSConcentration +=
221  conditions.SSConcentrations().back()
222  - *(conditions.SSConcentrations().end() - 2);
223  }
224  else
225  {
226  break;
227  }
228  }
229  dS = conditions.dV()
230  / (kdCalculator(currentSSConcentration) - 1.0)
231  / conditions.columnPoreVolume();
232  S += dS;
233  }
234 
235  RT = j * conditions.dV() / conditions.flowRate();
236  // Correction for the discreteness of integration.
237  if ((!backwardCompatibility) && (S > 1.0))
238  {
239  RT -= (S - 1.0) / dS * conditions.dV() / conditions.flowRate();
240  }
241  // Correction for Vp.
242  RT += conditions.columnPoreVolume() / conditions.flowRate();
243 
244  }
245 
246  RT += conditions.delayTime();
247  // Correction for V0.
248  RT += conditions.columnInterstitialVolume() / conditions.flowRate();
249  return RT;
250 }
251 }
252 
253 double calculateRT(const std::string &sequence,
254  const ChemicalBasis &chemBasis,
255  const ChromoConditions &conditions,
256  const int numInterpolationPoints,
257  const bool continueGradient,
258  const bool backwardCompatibility)
259  throw(BioLCCCException)
260 {
261  std::vector<ChemicalGroup> parsedSequence =
262  parseSequence(sequence, chemBasis);
263  return calculateRT(parsedSequence,
264  chemBasis,
265  conditions,
266  numInterpolationPoints,
267  continueGradient,
268  backwardCompatibility);
269 }
270 
271 double calculateKd(const std::string &sequence,
272  const double secondSolventConcentration,
273  const ChemicalBasis & chemBasis,
274  const double columnPoreSize,
275  const double columnRelativeStrength,
276  const double temperature)
277  throw(BioLCCCException)
278 {
279  return calculateKd(parseSequence(sequence, chemBasis),
280  secondSolventConcentration,
281  chemBasis,
282  columnPoreSize,
283  columnRelativeStrength,
284  temperature);
285 }
286 
287 double calculateAverageMass(const std::string &sequence,
288  const ChemicalBasis &chemBasis)
289  throw(BioLCCCException)
290 {
291  std::vector<ChemicalGroup> parsedSequence =
292  parseSequence(sequence, chemBasis);
293  double peptideAverageMass = 0;
294  for (std::vector<ChemicalGroup>::const_iterator i =
295  parsedSequence.begin();
296  i < parsedSequence.end();
297  i++)
298  {
299  peptideAverageMass += i -> averageMass();
300  }
301 
302  return peptideAverageMass;
303 }
304 
305 double calculateMonoisotopicMass(const std::string &sequence,
306  const ChemicalBasis &chemBasis)
307  throw(BioLCCCException)
308 {
309  std::vector<ChemicalGroup> parsedSequence =
310  parseSequence(sequence, chemBasis);
311  double peptideMonoisotopicMass = 0;
312  for (std::vector<ChemicalGroup>::const_iterator i =
313  parsedSequence.begin();
314  i < parsedSequence.end();
315  i++)
316  {
317  peptideMonoisotopicMass += i -> monoisotopicMass();
318  }
319 
320  return peptideMonoisotopicMass;
321 }
322 
323 }
double calculateAverageMass(const std::string &sequence, const ChemicalBasis &chemBasis)
Calculates the average (molar) mass of a peptide.
Definition: biolccc.cpp:287
An instance of ChemicalBasis contains a set of BioLCCC constants.
Definition: chemicalbasis.h:86
double calculateKd(const std::string &sequence, const double secondSolventConcentration, const ChemicalBasis &chemBasis, const double columnPoreSize=100.0, const double columnRelativeStrength=1.0, const double temperature=293.0)
Calculates the coefficient of distribution Kd for the given peptide.
Definition: biolccc.cpp:271
std::vector< ChemicalGroup > parseSequence(const std::string &source, const ChemicalBasis &chemBasis)
Parses a given peptide sequence.
Definition: parsing.cpp:134
Apart from classes, BioLCCC contains calculation methods and constants.
Definition: auxiliary.h:4
void fitSpline(const double *x, const double *y, const int n, double *y2)
Calculates the second derivatives of a function.
Definition: auxiliary.cpp:8
Base class for all BioLCCC exceptions. Can be used by itself.
double calculateSpline(const double *x, const double *y, const double *y2, const int n, const double x_in)
Calculates the value of a function using the cubic spline interpolation.
Definition: auxiliary.cpp:36
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
A ChromoConditions instance describes conditions of chromatography.
double calculateRT(const std::string &sequence, const ChemicalBasis &chemBasis, const ChromoConditions &conditions=standardChromoConditions, const int numInterpolationPoints=0, const bool continueGradient=true, const bool backwardCompatibility=false)
Calculates the retention time of a peptide.
Definition: biolccc.cpp:253
double calculateMonoisotopicMass(const std::string &sequence, const ChemicalBasis &chemBasis)
Calculates the monoisotopic mass of a peptide.
Definition: biolccc.cpp:305
double calculateKdRod(const std::vector< ChemicalGroup > &parsedSequence, const double secondSolventConcentration, const ChemicalBasis &chemBasis, const double columnPoreSize, const double columnRelativeStrength, const double temperature)
Calculates the coefficient of distribution of a polymer using the rod model.
Definition: rod_model.cpp:334