基于优化递归算法的分子量分解问题

基于优化递归算法的分子量分解问题摘要：本文讨论的问题是：在实验室拥有或不拥有计算机的情况下，如何将已知分子量x 的蛋白质分解成18种已知分子量的氨基酸的问题，并满足蛋白质含氮量在15%—17%的物理性质。

在实验室有计算机的情况下，本文首先考虑了穷举算法，该问题就等效为十八元一次方程求整数解的问题，表示为∑==181)](*)([i x i root i a每种氨基酸构成蛋白质数量root(i)的上限为])([i a x+1，取x 最大为1000时，有1610*8154.2种，运算量过大难以实现。

我们又考虑了递归算法,即从分子量最大的第18种氨基酸开始考虑，第18种氨基酸分子构成蛋白质的个数]])18([,0[)18(a xroot ∈，对root 进行分类取值，在取值已知的情况下再考虑第17种氨基酸，此时分子量减去第18种氨基酸总分子质量，得到新的分子量，将大大减少计算量。

抽象归纳为)(i root 的取值范围如下：)17,...,2,1]()())(*)(([)(0181=-≤≤∑+=i i a k root k a x i root i k该算法同样满足含氮量的约束条件。

我们借助C 语言编写程序并对程序进行优化，大大加快了运行速度。

在实验室无计算机的情况下，由氨基酸的结构通式得出，每个氨基酸都存在结构----CO CH NH ，分子量为56。

本文将18种氨基酸的分子量都减去56，得到18种氨基酸R 基的分子量。

先将任一分子量除以56取整，得到构成该蛋白质的氨基酸数目的范围（即]56[1xb ≤≤），再根据b 的取值进行分类讨论从而选择R 基，后根据蛋白质含氮量对结果进行验证，在人工的情况下，此方法相对于递归法大大地减少了计算量。

关键词：分子量分解问题递归优化法含氮量C语言结构通式1、问题重述生命蛋白质是由若干种氨基酸经不同的方式组合而成。

在实验中，为了分析某个生命蛋白质的分子组成，通常用质谱实验测定其分子量x (正整数)，然后将分子量x分解为n个已知分子量a[i](i=1,.......,n)氨基酸的和的形式。

某实验室所研究的问题中：n=18, x≤1000，a[i](i=1,.......,18)分别为57, 71, 87, 97, 99, 101, 103, 113, 114, 115, 128, 129, 131, 137, 147, 156, 163, 186。

要求针对该实验室拥有或不拥有计算机的情况，对如何分解分子量x作出解答，即针对任意一个分子量x具体给出由哪些a[i](i=1,.......,n)氨基酸组成。

2、问题分析题目中给出18种氨基酸分子量，要求解出已知分子量的蛋白质由哪几种氨基酸组成。

初步分析，这是十八元一次方程求整数解的问题，方程形式为∑==181)](*)([ixibia其中x为已知量，)18,...,2,1)((=iia为18种氨基酸分子量，未知数)18,...,2,1)((=iib为组成蛋白质的每种氨基酸的个数。

3、模型假设（1）假设氨基酸组合形成蛋白质时不脱水。

（2）假设蛋白质中肽键不存在多键或成环情况，两个氨基酸之间只形成一条肽键。

（3）假设每种氨基酸出现的概率都相同。

4、符号设定)(ia—第i种氨基酸的分子量；)(i root — 第i 种氨基酸的数量；)(_i number N — 第i 种氨基酸分子中氮原子的个数； x — 蛋白质分子量；)(i r — 第i 种氨基酸R 基的分子量； b — 蛋白质中所含氨基酸的总数；5、模型建立与求解5.1实验室有计算机时的模型建立5.1.1穷举法模型我们首先考虑用穷举算法求解方程，根据上文符号设定，该方程可以表示为∑==181)](*)([i x i root i a其中]])([,0[)(i a xi root ∈且)(i root 为整数。

利用穷举法，对每个)(i root 分别取0，1，…，])([i a x，找出符合上述方程的每个)(i root 的取值，所得的每组解的集合即为最终解。

穷举算法利用MATLAB 实现，程序代码见附录一。

用上述算法，我们发现，如果x 的取值太大，MATLAB 警告“busy ”，程序无法运行。

假如1000=x ，求出每一个)(i root 的最大值，如表一所示：表一 )(i root 最大值运算次数为1610*8154.2次，MATLAB 无法进行如此大数量级的运算，导致可求解的x 具有很大的局限性。

所以我们建立下文中的新模型，从而减少运算次数，提高运算速度，增强运算可行性。

5.1.2递归法模型为减少冗余次数，我们设计思路如下：从分子量最大的氨基酸开始考虑，初步确定第18种氨基酸分子个数后，对剩下17种氨基酸的确定，则不需要再用x 作为衡量标准，而可以将蛋白质中所含第18种氨基酸总分子质量减掉后，得到新的分子量，则剩下的17种氨基酸可以用新的分子量来确定。

按照此思路分析，)(i root 的新取值范围如下：)17,...,2,1]()())(*)(([)(0181=-≤≤∑+=i i a k root k a x i root i k利用这种模型，可以有效减少运算次数，使程序可以顺利运行。

5.1.3结果筛选我们考虑通过查阅资料了解到测定蛋白质的常用指标是含氮量，一般蛋白质的含氮量在15%~17%左右]1[。

所以我们设定限制条件%17*)(_*14%15*181x i number N x i ≤≤∑=利用此条件，可以把不符合蛋白质基本性质的结果筛选出去，剩下的解即为该问题的最终结果。

根据题目中给出的18种氨基酸的分子量，我们可以查找出题目中18种氨基酸的种类、分子式，从而得到每种氨基酸分子中氮原子个数)(_i number N ，如表二所示（题目中按照氨基酸两端成键给出分子量，故在写分子式时也考虑为两端成键）：表二 18种氨基酸信息一览表根据5.1.2中的递归模型以及5.1.3中的筛选条件，我们利用C语言编程实现，程序代码见附录二。

当261x时，用此程序计算结果如表三所示：=表三当261=x时利用递归法求解结果5.1.4优化递归算法在附录二的程序中，我们调用了递归函数求解，而层层调用递归函数往往会使预算速度变慢，不利于该方法的实际应用价值的发挥。

所以，我们对程序进行优化，不再调用递归函数，而取一个数列temp ，将上次运算后的结果储存在temp 的每一个格子中，这样可以避免多次调用递归函数导致的冗余运算，并将结果以txt 格式输出，不再结果框中显示结果，从而大大加快程序运行速度，优化后程序见附录3。

当1000=x 时，优化前运算时间为s 08.3，优化后运算时间不到s 1.0。

5.2实验室无计算机时的模型建立氨基酸结构通式为----CO CHR NH （两端成键，R 为任意元素或基），不难得出，每种氨基酸都存在基本结构----CO CH NH ，该基本结构分子量为56，故可以得到18种氨基酸R 基分子量。

先将任一分子量除以56取整，得到构成该蛋白质的氨基酸数目的范围（即]56[1xb ≤≤），再根据b 的取值进行分类讨论从而选择R 基，后根据蛋白质含氮量对结果进行验证，R 基分子量如表四所示：表四 18种氨基酸R 基分子量仍以x=261为例，因为456261=⎥⎦⎤⎢⎣⎡，所以下面分4种情况进行分类讨论：（1）b=1 ，表示蛋白质由一个氨基酸构成，存在一个R基等于261-56=205，r=130,故该情况不成立。

max（2）b=2，表示蛋白质由两个氨基酸构成，存在两个R基之和等于261-56*2=149，58+91=149，即261=114+147，含N量为26%，不符合要求。

故该情况也不成立。

（3）b=3,表示蛋白质由三个氨基酸构成，存在三个R基之和等于261-56*3=93：当261=57+57+147故存在三种情况：（4）b=4261-56*4=37，不存在任一情况使之成立。

6、模型分析在实验室有计算机的情况下，本文结合化学基本知识，采用递归算法，并用C语言编程。

与此同时，我们对C语言程序进行了优化，使得速度大大提升，方便用于实际计算与测量中。

在实验室无计算机的情况下，本文找出了氨基酸的结构通式，明确了氨基酸基本结构，将问题从氨基酸组合转化为R基组合，使计算量大大减小。

但当x较大时，人工方法难以避免较大的计算量。

7、模型推广在实际计算与测量中，蛋白质的物理性质和化学性质是不容忽视的。

本文中只利用蛋白质含氮量这一个条件对结果进行筛选，但得出的结果未必全部满足蛋白质的其他性质，故在实际生产生活中，如想要准确测定某种蛋白质的组成，还需要借助其他专业仪器，从多方面着手。

8、结论在实验室有计算机的情况下，本文建立递归算法模型，根据蛋白质含氮量在15%—17%的物理性质，对所得结果进行筛选。

并利用C语言编程，优化了程序，形成“优化递归算法”，大大加快运行速度。

在实验室无计算机的情况下，每个氨基酸都存在结构-CHNH，分子量为--CO-56。

所以本文得到18种氨基酸R基的分子量，再根据已知量x，按照所含氨基酸个数进行分类讨论。

此方法相对于递归法，大大地减少了计算量，方便在无计算机的情况下解出最终结果。

9、参考文献[1]Mark M.Meerscharert，《数学建模方法与分析》，北京：机械工业出版社，2005。

附录1：MATLAB程序（穷举法）S=zeros(1,18);for a=0:fix(x/57)for b=0:fix(x/71)for c=0:fix(x/87)for d=0:fix(x/97)for e=0:fix(x/99)for f=0:fix(x/101)for g=0:fix(x/103)for h=0:fix(x/113)for i=0:fix(x/114)for j=0:fix(x/115)for k=0:fix(x/128)for l=0:fix(x/129)for m=0:fix(x/131)for n=0:fix(x/137)for o=0:fix(x/147)for p=0:fix(x/156)for q=0:fix(x/163)for r=0:fix(x/186)if57*a+71*b+87*c+97*d+99*e+101*f+103*g+113*h+114*i+115*j+128*k+129*l+131*m+137*n+147*o+ 156*p+163*q+186*r==xS=[S;a b c d e f g h i j k l m n o p q r];endendendendendendendendendendendendendendendendendendendS附录2：C程序（递归法）#include <cstdio>#include <iostream>#include <fstream>using namespace std;//global settingint weight[18]={57, 71, 87, 97, 99, 101, 103, 113, 114, 115, 128, 129, 131, 137, 147, 156, 163, 186}; int root[18]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};int N_number[18]={1,1,1,1,1,1,1,1,2,1,2,2,1,1,3,1,4,1};int K=0;double root_number=0;double root_number2=0;ofstream fout("E://TestDir//cfiletest.txt");void Digui(int n,int x){int branch_number = x/weight[n];//printf("%d",branch_number);if(x<0){return;}if(n == 0){int i = branch_number;for(i = branch_number;i >= 0;i--){root[n] = i;if( x-i*weight[0] == 0){double N_content=0;int k=0;for(k=0;k<18;k++){N_content=N_content+N_number[k]*root[k]*14;}//printf("K:%d\n",K);//printf("K*0.15:%f\n",K*0.15);//printf("%f\n",N_content);//printf("K*0.17:%f\n",K*0.17);if( N_content>K*0.15 && N_content<K*0.17){root_number = root_number + 1;int j=0;printf("root:");fout<<"root"<<root_number<<":";for(j=0;j<18;j++){printf("%d ",root[j]);fout<<root[j]<<" ";}printf("\n");fout<<endl;}root_number2=root_number2+1;}}return;}else{int i = branch_number;for(i = branch_number;i >= 0;i--){root[n] = i;//printf("x: %d\n",x-i*weight[n]);//printf("n-1:%d\n",n-1);Digui(n-1,x-i*weight[n]);}}return;}int main(){int x;printf("Please input K:");scanf("%d",&K);x=K;int i=0;printf("weight:");for(i=0;i<18;i++){printf("%d ",weight[i]);}printf("\n");Digui(17,x);printf("total:");printf("%f \n",root_number);printf("unstrict total:");printf("%f \n",root_number2);return 0;}附录3：C程序（优化递归算法）#include <stdio.h>//global settingint weight[18]={57, 71, 87, 97, 99, 101, 103, 113, 114, 115, 128, 129, 131, 137, 147, 156, 163, 186}; int root[18]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};int x=0;double root_number=0;int temp[18];int main(){printf("Please input x:");scanf("%d",&x);int i;int j[18];for(j[17]=x/weight[17];j[17]>=0;j[17]--){temp[17]=x-weight[17]*j[17];for(j[16]=temp[17]/weight[16];j[16]>=0;j[16]--){temp[16]=temp[17]-weight[16]*j[16];for(j[15]=temp[16]/weight[15];j[15]>=0;j[15]--){temp[15]=temp[16]-weight[15]*j[15];for(j[14]=temp[15]/weight[14];j[14]>=0;j[14]--){temp[14]=temp[15]-weight[14]*j[14];for(j[13]=temp[14]/weight[13];j[13]>=0;j[13]--){temp[13]=temp[14]-weight[13]*j[13];for(j[12]=temp[13]/weight[12];j[12]>=0;j[12]--){temp[12]=temp[13]-weight[12]*j[12];for(j[11]=temp[12]/weight[11];j[11]>=0;j[11]--){temp[11]=temp[12]-weight[11]*j[11];for(j[10]=temp[11]/weight[10];j[10]>=0;j[10]--){temp[10]=temp[11]-weight[10]*j[10];for(j[9]=temp[10]/weight[9];j[9]>=0;j[9]--){temp[9]=temp[10]-weight[9]*j[9];for(j[8]=temp[9]/weight[8];j[8]>=0;j[8]--){temp[8]=temp[9]-weight[8]*j[8];for(j[7]=temp[8]/weight[7];j[7]>=0;j[7]--){temp[7]=temp[8]-weight[7]*j[7];for(j[6]=temp[7]/weight[6];j[6]>=0;j[6]--){temp[6]=temp[7]-weight[6]*j[6];for(j[5]=temp[6]/weight[5];j[5]>=0;j[5]--){temp[5]=temp[6]-weight[5]*j[5];for(j[4]=temp[5]/weight[4];j[4]>=0;j[4]--){temp[4]=temp[5]-weight[4]*j[4];for(j[3]=temp[4]/weight[3];j[3]>=0;j[3]--){temp[3]=temp[4]-weight[3]*j[3];for(j[2]=temp[3]/weight[2];j[2]>=0;j[2]--){temp[2]=temp[3]-weight[2]*j[2];for(j[1]=temp[2]/weight[1];j[1]>=0;j[1]--){temp[1]=temp[2]-weight[1]*j[1];for(j[0]=temp[1]/weight[0];j[0]>=0;j[0]--){if(temp[1]-weight[0]*j[0]==0){root_number++;}}}}}}}}}}}}}}}}}}}printf("%f\n",root_number);return 0;}。