text:
# Generalized-Ensemble Simulations of the Human Parathyroid Hormone Fragment PTH(1-34)

## I Introduction

The successful deciphering of the human genome has aggravated an old challenge in protein science: for most of the resolved protein sequences we do not know the corresponding structures and functions. Computer experiments offer one way to evaluate the sequence-structure relationship but are extremely difficult for realistic protein models where interactions among all atoms are taken into account. The complex form of the intramolecular forces and of the interaction with the solvent, containing both repulsive and attractive terms, leads to a very rough energy landscape with a huge number of local minima. These minima are separated by energy barriers that are much higher than the typical thermal energy of a protein (of order $`k_BT`$) at room temperature. Hence, simple canonical Monte Carlo or molecular dynamics simulations will get trapped in a local minimum and often not thermalize within a finite amount of available CPU time. While this multiple minima problem does not necessarly inhibits molecular dynamics or Monte Carlo studies of peptides and proteins Calfisch1 ; Kol , it restricts calculation of accurate thermodynamic averages to small peptides Chow .

A number of novel simulation techniques have been proposed for overcoming this multiple-minima problem (for a review, see Ref. HO98g, ). Important recent examples can be found in Ref. Brooks2002, ; Elber, . Another example is parallel tempering, also known as replica exchange method and introduced to protein science in Ref. H97f, , that has become increasingly popular over the last few years Angel ; LHH ; Calfisch2 . Parallel tempering is only one example of a class of new and sophisticated algorithms commonly summarized as generalized-ensemble methods MyReview . In this article, we are concerned with another generalized-ensemble technique, multicanonical sampling MU , that was first applied to the protein-folding problem in Ref. HO1, . Its usefullnes for calculating reliable thermodynamic quantities at low temperatures has been established for small peptides of up to $`20`$ residues HO1 ; HO98c ; HA02 . However, stable domains in proteins consists usually of 40-200 amino acids. Hence, it is necessary to evaluate the effectiveness of this approach for larger molecules than the peptides investigated so far. For this purpose, we have performed generalized-ensemble simulations of the peptide fragment PTH(1-34) corresponding to residues 1-34 of human parathyroid hormone Klaus ; Crystal ; MABFR .

The 84-amino acid human parathyroid hormone is involved in the regulation of the calcium level in blood and influences bone formation Potts . The NH2-terminal 34 residues of the hormone, further on refered to as PTH(1-34), are sufficient for the biological activities of this hormone suggesting medical and pharmaceutical applications of the peptide Brommage99 . The crystal structure of the peptide has been resolved at 0.9 Å and resembles a slightly bend long helix (PDB code 1ET1) Crystal . NMR studies of the peptide in solution under near physiological conditions (PDB code 1ZWA) MABFR  and in 20 % trifluorethanol solution (PDB code 1HPY) MABFR  rather indicate an ensemble of structures that have in common two helices separated by a disordered region.

In the present article, we try to overcome the problems of previous simulated annealing simulations of PTH(1-34) Okamoto , that did not allow a detailed structure evaluation, by using multicanonical sampling MU , one of the most prominent generalized-ensemble techniques. An all-atom representation of the molecule is employed and the intramolecular interactions are described by the ECEPP/3 force field EC3 . The protein-solvent interactions are approximated by the solvent accessible surface term of Ooi et al.OONS  Quantities such as the average helicity, number of contacts, average energy, and specific heat are calculated. Our results demonstrate the feasibility of generalized-ensemble simulations for large molecules such such as PTH(1-34). In addition, they indicate that with the advent of these and other modern search techniques, structure prediction of proteins is limited more by current energy functions (and especially solvent approximations) than by the simulation algorithms.

## II Methods

Our research into the thermodynamics of PTH(1-34) is based on a detailed, all-atom representation of that peptide. The interactions between the atoms are described by a standard force field, ECEPP/3 EC3  (as implemented in the program package SMMP SMMP ), and are given by the sum of the electrostatic term $`E_C`$, the Lennard-Jones energy $`E_{LJ}`$, hydrogen-bond term $`E_{HB}`$ for all pairs of atoms in the peptide together with the torsion term $`E_{tor}`$ for all torsion angles:

$`E_{ECEPP/2}`$ $`=`$ $`E_C+E_{LJ}+E_{HB}+E_{tor},`$ (1)
$`E_C`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}{\displaystyle \frac{332q_iq_j}{ϵr_{ij}}},`$ (2)
$`E_{LJ}`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}\left({\displaystyle \frac{A_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{B_{ij}}{r_{ij}^6}}\right),`$ (3)
$`E_{HB}`$ $`=`$ $`{\displaystyle \underset{(i,j)}{}}\left({\displaystyle \frac{C_{ij}}{r_{ij}^{12}}}{\displaystyle \frac{D_{ij}}{r_{ij}^{10}}}\right),`$ (4)
$`E_{tor}`$ $`=`$ $`{\displaystyle \underset{l}{}}U_l\left(1\pm \mathrm{cos}(n_l\chi _l)\right).`$ (5)
Here, $`r_{ij}`$ (in Å) is the distance between the atoms $`i`$ and $`j`$, $`\chi _l`$ is the torsion angle for the chemical bond $`l`$ and $`n_l`$ characterizes its symmetry. The charges and force field parameters $`q_i,A_{ij},B_{ij},C_{ij},D_{ij},U_l`$ were calculated from crystal structures of amino acids using semi-empirical methods. The dielectricity constant is set to $`\epsilon =2`$, its common value in ECEPP calculations. Since the bond lengths are fixed in ECEPP, the backbone torsion angles $`\varphi ,\psi ,\omega `$ and the side chain torsion angles $`\chi `$ are the true degrees of freedom. Hence, in a Monte Carlo (MC) sweep single angles are updated sequentially by the Metropolis algorithm. The protein-water interactions are approximated by a solvent-accessible surface term

$$E_{solv}=\underset{i}{}\sigma _iA_i,$$
(6)
where $`A_i`$ is the solvent accessible surface area of the $`ith`$ atom in the present configuration, and $`\sigma _i`$ the solvation parameter for the atom $`i`$. We choose the solvation parameter set OONS of Ref. OONS  that is often used together with the ECEPP force field. The potential energy of the solvated molecule is then given by

$$E_{tot}=E_{ECEPP/3}+E_{solv}$$
(7)
In such a detailed protein model, the various competing interactions lead model to an energy landscape with a multitude of local minima separated by high-energy barriers. Canonical Monte Carlo or molecular dynamics simulations will likely get trapped in one of these minima and not thermalize within the available CPU time. Only recently, with the introduction of new and sophisticated algorithms such as generalized-ensemble techniques MyReview  was it possible to alleviate this problem in protein simulations HO1 . For simulating PTH(1-34) we have chosen one most commonly used generalized-ensemble technique, multicanonical sampling MU .

The multicanonical algorithm MU  assigns a weight $`w_{mu}(E)1/n(E)`$ to conformations with energy $`E`$. Here, $`n(E)`$ is the density of states. A simulation with this weight leads to a uniform distribution of energy:

$$P_{mu}(E)n(E)w_{mu}(E)=\mathrm{const}.$$
(8)
Thus, the simulation generates a 1D random walk in the energy space, allowing itself to escape from any local minimum. Since a large range of energies is sampled, re-weighting FS  allows one to calculate thermodynamic quantities over a wide range of temperatures $`T`$ by

$$<𝒜>_T=\frac{𝑑x𝒜(x)w^1(E(x))e^{\beta E(x)}}{𝑑xw^1(E(x))e^{\beta E(x)}},$$
(9)
where $`x`$ stands for configurations, $`E(x)`$ for its total potential energy $`E(x)=E_{ECEPP/3}(x)+E_{solv}(x)`$ and $`\beta `$ for the inverse temperature, $`\beta =1/k_BT`$.

Unlike in constant temperature simulations the weights are not a priori known in multicanonical simulations. In fact, knowledge of the exact weights is equivalent to obtaining the density of states $`n(E)`$, i.e., solving the system. However, for a numerical simulations estimators are sufficient as long as these do not deviate not too much from $`n^1(E)`$. This is because the same weights that are used for the simulation appear also in the re-weighting procedure of Eq. 9. In the present study, we calculate these estimators from a preliminary simulated annealing run of 80,000 MC-sweeps using the method described in Ref. H97, . All thermodynamic quantities are then estimated from one production run of $`1,000,000`$ sweeps, starting from a random initial configuration and after discarding $`10,000`$ sweeps for thermalization. We store in every fifth sweep for further analysis various physical quantities and the dihedral angles of the current configuration. Our error bars are estimated by dividing this time series of data into 8 bins of each $`125,000`$ sweeps.

## III Results and Discussion

We start our analysis by calculating thermodynamic averages of the intramolecular energy $`<E_{ECEPP/3}>(T)`$ and the solvation energy $`<E_{SOLV}>(T)`$. Both quantities and the resulting total energy $`<E_{TOT}>(T)=<E_{ECEPP/3}+E_{SOLV}>(T)`$ are displayed as function of temperature in Fig. 1. The thermal behavior of the peptide is characterized by a competition between intramolecular and solvation energy. While $`<E_{ECEPP/3}>(T)`$ decreases with decreasing temperature, $`<E_{SOLV}>(T)`$ increases toward lower temperatures. The interplay of both terms leads to two temperature regimes that are separated by a steep decrease in the total energy $`<E_{TOT}>(T)`$. The corresponding pronounced peak in the specific heat per residue

$$C(T)=\beta ^2\left(<E_{TOT}^2>(T)<E_{TOT}>^2(T)\right)/34,$$
(10)
displayed in the inlet, marks the transition temperature at $`T_c=560\pm 10`$ K.

The structural changes associated with this transition can be deduced from Fig. 2 where we display the average helicity $`<n_H>(T)`$ as a function of temperature. Here, we have defined $`n_H`$ as the number of residues whose pair of backbone dihedral angles $`(\varphi ,\psi )`$ takes values in the range: $`(70^{}\pm 30^{},37^{}\pm 30^{})`$. We see from the plot of both quantities that the high-temperature (high-energy) region is characterized by configurations with vanishing helicity ($`10`$%) while at low temperatures (and, correspondingly, low energies) configurations dominate that are almost completely helical ($`90`$% of the residues are part of a helix). The pronounced peak at $`T_c`$ in the susceptibility (per residue)

$$\chi (T)=\left(<n_H^2><n_H>^2\right)/34,$$
(11)
shown in the inlet, is further proof for the sharp transition between low-energy helical states and high-energy disordered coil states. Associated with this helix-coil transition is also a decrease in the solvent accessible volume $`<V>(T)`$ as calculated by the double cubic lattice method ELASS  (data not shown). Above $`T_c`$, the average volume is $`<V>10000`$ Å3, while below $`T_c`$ the volume is reduced to $`<V>9000`$ Å3.

The modest decrease in $`<V>`$ together with the large value of the helicity $`<n_H>(T)`$ indicate that below $`T_c`$ a single elongated helix is the dominant structure for PTH(1-34). Indeed, we find in our simulation as lowest-energy state an elongated helix with 31 residues part of the helix. This structure, shown in Fig. 3b, has not only the lowest total energy ($`E_{TOT}=277.8`$ kcal/mol), but also the lowest intramolecular energy: $`E_{ECEPP/3}=136.5`$ kcal/mol. It is very similar to the crystal structure of PTH(1-34) (PDB code 1ET1, displayed in Fig. 3a) where also 31 residues are part of an $`\alpha `$-helix and whose energy after regularization with the program FANTOM FANTOM  is $`E_{TOT}=277.9`$ kcal/mol ($`E_{ECEPP/3}=187.0`$ kcal/mol). While our numerically determined structure has a slightly larger solvent-accessible surface area ($`A=3860`$ Å2) than the crystal structure ($`A=3410`$ Å2), it has 95% of all native contacts formed, i.e. 95 % of the contacts between residues found in the crystal structure exist also in our lowest-energy configuration. Here, we consider two residues in contact if the distance between their Cα-atoms is less than 8 Å, and the two residues are neither neighbor nor next-nearest neighbor in the peptide chain. Given that almost all native contacts are formed in the lowest-energy structure, it is not surprising that the root-mean-square deviation (rmsd) between the two structures is only $`0.8`$ Å for backbone atoms ($`2.3`$ Å when all heavy atoms are taken into account). For comparison, the crystal structure of PTH(1-34) (1ET1) itself was solved at $`0.9`$ Å resolution Crystal . We remark that recent structure determinations of the similar sized villin headpiece subdomain HP-36, a 36-residue peptide, by Energy Landscape Paving HW  and parallel tempering LHH  were restricted to an accuracy of $`6`$ Å. We believe that the higher accuracy of our PTH(1-34) results does not indicate any advantage of multicanonical sampling over the above methods but is rather due to the simpler geometry of PTH(1-34).

At $`T=300`$ K, our lowest energy structure appears with a frequency of $`(99\pm 0.5)`$%, i.e. almost all observed configurations resemble the crystal structure. The predominance of this structure is also reflected in the well-developed funnel in Fig. 4 where we display the projection of the free-energy landscape at $`T=300`$ K on the number of native contacts. The free energy decreases rapidly with increasing number of native contacts. No indications for competing local minima that could act as traps are observed indicating a rather smooth funnel. On the other hand, at the transition temperature $`T_c=560`$ K, the free energy landscape (displayed in the inlet) is flat and configurations with small number of native contacts coexist with such that have many native contacts.

However, while the crystal structure is in our simulation the dominant configuration at $`T=300`$ K, it differs from the set of NMR-structures found at room temperature. In near-physiological solution, one observes instead two helices separated by a disordered and flexible region. We show in Fig. 3c as an example one of the resolved solution configurations (from 1HPY) MABFR . The N-terminal helix ranges from Glu4 to His9 and the C-terminal helix from Ser17 to Gln29. Addition of trifluorethanol reduces hydrophobic interactions and increases the length of these helices MABFR . Hence, our simulation of PTH(1-34) does not reproduce the experimental results for that peptide in solution albeit protein-solvent interactions are considered in our energy function by an approximate term. Instead, our simulation favors the crystal structure of the peptide that is observed in membrane and hydrophobic environments.

In order to understand in greater detail the relation between our simulation results and the NMR experiments, we plot in Fig. 5a for each residue the free energy difference $`\mathrm{\Delta }G_i`$ at $`T=300`$ k between configurations with residue $`i`$ part of an $`\alpha `$-helix and such where that residue is not part of an $`\alpha `$-helix. The free-energy differences are largest for residues Asn16 -Lys27, and it is for these residues that first helix formation is observed. A second cluster of residues that have large free-energy differences are observed between Ile5 and Asn10. Both regions are separated by residues Leu11 \- Leu15 that have smaller free energy differences. The observed free-energy differences are strongly correlated with differences in the (potential) energy $`\mathrm{\Delta }E_{TOT}`$ that together with its two components $`E_{ECEPP/3}`$ and $`E_{SOLV}`$ are plotted in Fig. 5b. Note, that the variations in the energy differences result from the $`E_{ECEPP/3}`$ part, i.e. from the intramolecular interactions. The corresponding solvent energies favor in general residues that are not in a helical state but vary little with the residues. In addition, their magnitude is so small that it is difficult to distinguish in the figure between $`\mathrm{\Delta }E_{TOT}`$ and $`\mathrm{\Delta }E_{ECEPP/3}`$.

The position of the two helices in the solvent structure of Fig. 3c corresponds to the regions where in our simulation the measured free-energy differences and potential energy differences are large. In the NMR structures, the C-terminal helix is more stable than the N-terminal helix. Similarly, we find larger absolute values of $`\mathrm{\Delta }G_i14`$ kcal/mol ($`\mathrm{\Delta }E_{TOT}39`$ kcal/mol) for residues Asn16 to Lys27 (with the maximal values at Arg20: $`\mathrm{\Delta }G=19.4`$ kcal/mol and $`\mathrm{\Delta }E_{TOT}=55.1`$ kcal/mol) compared with $`\mathrm{\Delta }G_i12`$ kcal/mol ($`\mathrm{\Delta }E_{TOT}36`$ kcal/mol) for residues Ile5 to Asn10. On the other hand, the flexible region connecting the two helices in the NMR structure corresponds to a region of residues that have with $`\mathrm{\Delta }G_i8`$ kcal/mol and $`\mathrm{\Delta }E_{TOT}25`$ kcal/mol considerably smaller free (potential) energy differences. The free-energy differences are smallest for Leu11 and Gly12: $`\mathrm{\Delta }G_i6`$ kcal/mol and and $`\mathrm{\Delta }E_{TOT}17`$ kcal/mol. The later result is not surprising giving the inherent flexibility of glycine (which, however, is part of the helix in our lowest-energy configuration).

The observed variations in the free and potential energy differences suggest that we may find at higher temperatures configurations similar to the NMR structures. This is because the helix will be de-stabilized with increasing temperature, and more easily for residues Leu11 to Leu15 than in the regions that corresponds to the two terminal helices. In order to test this conjecture, we show in Fig. 6 two quantities as function of temperature. One is the frequency of configurations that have a continuous helix stretching at least between Ile5 and Lys27 and are therefore similar to the crystal structure of PTH(1-34). The second quantity is the frequency of configurations that have helices stretching at least between Ile5 and His9 and between Ser17 and Gln29, but are separated in-between by a non-helical flexible region. Hence, the later quantity measures the frequency of configurations that are similar to the NMR structure. A typical example of this group of configurations is displayed in Fig. 3d. While this conformation is also a local minimum, its total energy $`E_{TOT}=201.4`$ kcal/mol and intramolecular energy $`E_{ECEPP/3}=48.6`$ kcal/mol are much higher than the corresponding values for the lowest-energy configuration of Fig. 3b: $`E_{TOT}=277.8`$ kcal/mol and $`E_{ECEPP/3}=135.5`$ kcal/mol. On the other hand, its solvation energy $`E_{SOLV}=153.0`$ kcal/mol is lower than that of Fig. 3b ($`E_{SOLV}=141.3`$ kcal/mol) but the differences are smaller than the one in the ECEPP/3 term.

Both kind of configurations appear at the helix-coil transition temperature $`T_c=560(10)`$ K. Due to their higher entropy, configurations that resemble the NMR structures are slightly more common for temperatures in the range $`520`$ K $`T560`$ K than the ones that are similar to the crystal structure leading to a positive free energy difference $`\mathrm{\Delta }G`$ that is displayed in the inlet of Fig. 6. At $`T=520`$ K, 43(4)% of all configurations are similar to the NMR structures and 40(5)% resemble more the crystal structure. Below $`T=520`$ K, the intramolecular energy that favors an extended single helix wins over the higher entropy of states with two separate helices. The resulting negative free energy difference $`\mathrm{\Delta }G`$ leads to a decrease in the frequency of NMR-like structures, and their contribution is less than 1% at room temperature.

We conjecture that the above relations hold also in nature. The elongated helix of the crystal structure is favored by the intramolecular energies and is the ground state in potential energy. Thermal fluctuations lead to the more flexible configurations with two helices that are observed for the soluted peptide. In a (hydrophobic) membrane environment or when binding to a receptor reduces the entropy of the molecule, PTH(1-34) stays in the 1-helix state that also seems to lead to increase its biological activity Crystal .

However, the energy function in our simulation overstabilizes the extended $`\alpha `$-helix of the PTH(1-34) ground state (the crystal structure). Hence, structures that are found in solution with NMR experiments appear in our simulation with significant frequency only at temperatures more than $`200`$ K above room temperature. Similarly, we find with $`T_c=560(10)`$ K a helix-coil transition temperature that is much higher than physiological relevant temperature range. These results clearly point out the limitations of our model. Since our data are consistent with what one would expect in a hydrophobic environment one can conjecture that these limitations of the energy function are mainly due to our solvent approximation. The small variation in the magnitude of the solvent energy in Figs. 1 and 5 (when compared with the ECEPP/3 term) suggests that the OONS term underestimates the protein-solvent interaction. This maybe due to the fact that our implicit solvent model does not account for the hydrogen bonds between the polypeptide and water molecules that in water compete with the characteristic intramolecular hydrogen bonding in an $`\alpha `$-helix. As a consequence, $`\alpha `$-helices are overstabilized, and our energy function rather models the peptide in a hydrophobic environment than in physiological solution. Another problem is that the solvation energy term of Eq. 6 describes actually a free energy and therefore should change with temperature. This effect is neglected in the OONS approximation. Taking such a temperature dependence into account and use of more sophisticated implicit solvent models would likely improve our results. However, we can also not exclude the possibility that the deviation from the NMR results is not duee to the solvent term bit that our force field, the ECEPP/3 term that describes the intramolecular interactions, biases toward helical conformations.

It follows that other potential energy functions and implicit solvent models have to be chosen for a simulation of PTH(1-34) in solution. However, the failure of our energy function to model correctly the solvated molecule demonstrates also the advantages of our approach. Unlike the earlier simulated annealing simulations of Ref. Okamoto,  multicanonical sampling allows one to calculate accurate estimates of the frequency of states and other physical quantities at room temperature. In that way, we have been able to unveil in the present paper the limitations of our energy function. At the same time, our results allow us to understand the behavior of our peptide in a hydrophobic environment (since this is what our energy function models) and to reproduce the crystal structure with high accuracy. Hence, we have shown that multicanonical simulations are not only well-suited for simulations of small peptides but also for numerically more challenging molecules such as the 34-residue peptide PTH(1-34).

## IV Conclusion

In summary, we have performed all-atom simulations of PTH(1-34), the biologically active peptide-fragment 1-34 of the human parathyroid hormone. Protein-water interactions are approximated by a solvent-accessible surface term using the OONS parameter setOONS . Our results rely on multicanonical sampling and demonstrate that this technique allows one to overcome the multiple minima problem in simulations of this molecule. Unlike earlier numerical studies Okamoto  that relied on simulated annealing we find a lowest-energy configurations that has a rmsd of only 0.8 Å to the crystal structure. While previous applications of multicanonical sampling were restricted to small peptides with less than 20 residues, our results (albeit for a molecule with a rather simple structure) establish that the generalized-ensemble approach is also a useful tool for investigation of much larger polypeptides with $`3040`$ residues. Applications to even larger molecules seem to be restricted less by the sampling technique but by the accuracy of the energy function.

Acknowledgments:
Support by a research grant (CHE-9981874) of the National Science Foundation is gratefully acknowledged. Part of the described research was done while I was visiting the John von Neumann Institute for Computing (NIC), located in the FZ Jülich. I like to thank Prof. Grassberger and the institute for kind hospitality and for bearing with my working style.

Figure Captions:

Average intramolecular energy $`<E_{ECEPP/3}>`$, solvent energy $`<E_{SOLV}>`$ and total energy $`<E_{TOT}>=<E_{ECEPP/3}+E_{SOLV}>`$ of PTH(1-34) as a function of temperature $`T`$. The inlet displays the specific heat $`C(T)`$ as a function of temperature $`T`$. The data are calculated from a multicanonical simulation of 1,000,000 sweeps using a solvent accessible surface term to approximate protein-water interactions.

Average helicity $`<n_H>`$ as a function of temperature $`T`$. The inlet displays the susceptibility $`\chi (T)`$ as a function of temperature $`T`$. All data points are calculated from a multicanonical simulation of 1,000,000 sweeps using a solvent accessible surface term to approximate protein-water interactions.

(a) The crystal structure of PTH(1-34) (PDB code 1ET1); (b) The lowest-energy conformation of PTH(1-34) as determined from a multicanonical simulation with a solvent accessible surface term; (c) Solution structure of PTH(1-34) as determined by NMR (PDB-code 1HPY). (d) One of the configurations resembling the solution structure. Such configurations appear in our simulation with significant frequency only for temperatures closely to but below the helix-coil transition temperature $`T_c=560`$ K.

Projection of the free-energy landscape ($`\mathrm{\Delta }G`$) of PTH(1-34) at $`T=300`$ K on the number of native contacts $`n_{NC}`$. The same quantity is plotted for the $`T=560`$ K (the critical temperature) in the inlet.

(a) Free energy difference $`\mathrm{\Delta }G_i`$ between configurations with residue $`i`$ part of an $`\alpha `$-helix and such where residue $`i`$ is not in a helical state. (b) The corresponding total (potential) energy differences $`\mathrm{\Delta }E_{TOT}`$ ($`\mathrm{}`$), intramolecular energy differences $`\mathrm{\Delta }E_{ECEPP/3}`$ ($``$) and solvent energy differences $`\mathrm{\Delta }E_{SOLV}`$($`x`$).

Frequency of PTH(1-34) configurations ($`1H`$) that resemble the crystal structure and frequency of states ($`2H`$) that are similar to the solution structures of the peptide. The free-energy difference $`\mathrm{\Delta }G`$ between both sets of configurations is displayed in the inlet. All quantities are shown as function of temperature $`T`$. All our data points are calculated from a multicanonical simulation of 1,000,000 sweeps using a solvent accessible surface term to approximate protein-water interactions.


format: text

text:
# Do Pareto-Zipf and Gibrat laws hold true? An analysis with European Firms

## 1 Introduction

Pareto  is generally credited with the discovery, more than a century ago, that the distribution of personal income obeys a power-law in high-income range111See  for modern and high-quality personal-income data in Japan.. Firm size also has a skew distribution , and quite often obeys a power-law in the upper tail of the distribution. In terms of cumulative distribution $`P_>(x)`$ for firm size $`x`$, this states that

$$P_>(x)x^\mu ,$$
(1)
for large $`x`$, with $`\mu `$ being a parameter called Pareto index. The special case $`\mu =1`$ is often referred to as Zipf’s law . In this paper we call it Pareto-Zipf law, the fact that firm size has a power-law distribution asymptotically for large firms.

Even if the range for which eq. (1) is valid is a few percent in the upper tail of the distribution, it is often observed that such a small fraction of firms occupies a large amount of total sum of firm sizes. This means that a small idiosyncratic shock can make a considerable macro-economic impact. It is, therefore, quite important to ask what is the underlying dynamics that governs the growth of those large firms.

Let a firm’s size be $`x_1`$ at a time and $`x_2`$ at a later time. Growth rate $`R`$ is defined as the ratio $`R=x_2/x_1`$. Law of proportionate effect  (see also ) is a postulate that the growth rate of a firm is independent of the firm’s attained size, i.e.

$$P(R|x_1)\text{ is independent of }x_1,$$
(2)
where $`P(R|x)`$ is the probability distribution of growth rate conditional on the initial size $`x_1`$. In this paper we call this assumption as Gibrat’s law222Another interesting and related quantity is flow, e.g. profits, rather than stock. See  for growth of individual personal-income and  for firms tax-income growth, and validity of Gibrat’s law..

These two laws have been extensively studied in industrial organization and related stochastic models  (see  for review). Recent study in econophysics  introduced some notions and concepts of statistical physics into economics (see ). Present status related to firm-size growth may be summarized as follows. Firm size distribution is approximately log-normal with deviation from it in the upper tail of the distribution (e.g.  for recent data). On the other hand, Gibrat’s law breaks down in the sense that the fluctuations of growth rate scale as a power-law with firm size; smaller firms can possibly have larger fluctuations (e.g. ). However, little attention has been paid to the regime of firm size where power-law is dominant rather than log-normality, and to the validity of Gibrat’s law in that regime. More importantly, any kinematic relationship between Pareto-Zipf and Gibrat laws has not been understood explicitly, although there have been a lot of works on stochastic dynamics since Gibrat. This issue is exactly what the present paper addresses.

For our purpose it is crucial to employ exhaustive lists of large firms. Our dataset for European countries is exhaustive in the sense that each list includes all the active firms in each country whose sizes exceed a certain threshold of observation. We show that both of the Pareto-Zipf law and Gibrat’s law do hold for those large firms. As our main result, we prove that Pareto-Zipf law implies Gibrat’s law and vice versa under detailed balance. By showing that the condition of detailed balance also holds in our empirical data, we can show the equivalence of Pareto-Zipf law and Gibrat’s law as a kinematic principle in firms growth, irrespective of the underlying dynamics. Thereby, we conjecture that Gibrat’s law does hold in the regime of Pareto-Zipf for large firms, but does not for smaller firms. Thus our result is not contradictory to the breakdown of Gibrat’s law in previous study, most notably to the recent work by Stanley’s group . Furthermore, in the process of our proof, we also show that the distribution of growth rate possesses a non-trivial relation between the positive side ($`R>1`$) of the distribution and the negative side ($`R<1`$), through the value of Pareto index $`\mu `$, which is confirmed empirically.

In section 2, we give a brief review of the study on Gibrat’s law and firm size distribution in economics. In section 3, we describe the nature of our database of firms with large size in European countries. In section 4, using exhaustive lists of large firms in the dataset, we show that Gibrat’s law holds in the power-law regime for which the firm size distribution obeys Pareto-Zipf law. In addition, we uncover that temporal change of individual firm’s size in successive years satisfies what we call time-reversal symmetry, or detailed balance. In section 5, we prove that the two empirical laws of Gibrat and Pareto-Zipf are equivalent under the condition of detailed balance. We summarize our results in section 6.

## 2 Gibrat and Pareto-Zipf Laws in Economics

Industrial organization literature has long been focused on two empirical facts  (see  for review):
(i) skew distribution of firms size
(ii) validity or invalidity of Gibrat’s law for firm growth

Gibrat formulated the law of proportionate effect for growth rate to explain the empirically observed distribution of firms. The law of proportionate effect states that the expected increment to a firm’s size in each period is proportional to the current size of the firm. Let $`x_t`$ and $`x_{t\mathrm{\Delta }t}`$ be, respectively, the size of a firm at time $`t`$ and $`t\mathrm{\Delta }t`$, and $`ϵ_t`$ denote the proportionate rate of growth. The the postulate is expressed as

$$x_tx_{t\mathrm{\Delta }t}=ϵ_tx_{t\mathrm{\Delta }t}.$$

Gibrat assumed (a) that $`ϵ_t`$ is independent of $`x_t`$ (Gibrat’s law), (b) that $`ϵ_t`$ has no temporal correlation, and (c) that there is no interaction between firms. Then, after a sufficiently long time $`t\mathrm{\Delta }t`$, since

$$x_t=x_0(1+ϵ_1)(1+ϵ_2)\mathrm{}(1+ϵ_t),$$

$`\mathrm{log}x_t`$ follows a random walk. Assuming that $`ϵ_t`$ is small, one has

$$\mathrm{log}x_t=\mathrm{log}x_0+ϵ_1+ϵ_2+\mathrm{}+ϵ_t.$$

Gibrat’s model has two consequences concerning the above points (i) and (ii). Since the growth rate defined by $`R_tx_t/x_0`$ has its logarithm as the sum of independent variables $`ϵ_t`$, the growth rate is log-normally distributed. In addition, assuming that all the firms have approximately the same starting time and size, the distribution of firms size is also also log-normal with mean and variance given by $`mt`$ and $`\sigma ^2t`$, respectively, where $`m`$ is the mean of $`ϵ_t`$ and $`\sigma ^2`$ is the variance of $`ϵ_t`$.

The assumptions (a)–(c) in Gibrat’s model are in disagreement with empirical evidence. Among others, the Gibrat’s law (a) is incompatible with the fact that the fluctuations of growth rate measured by standard deviation decreases as firm size increases . Especially, the recent work  by Stanley’s group showed that the distribution of the logarithm of growth rates, for each class of firms with approximately the same size, displays an exponential form (Laplace distribution) rather than log-normal. They also show that the fluctuations in the growth rates characterized by the standard deviation $`\sigma (x)`$ of the distribution decreases for larger size of firms as a power-law, $`\sigma (x)x^\beta `$, with the exponent $`\beta `$ is less than a half. The latter point suggests a new viewpoint about the interplay of different parts of a firm, an industrial sector, or an organization .

In contrast to the standard deviation, the measure by mean growth rate has been disputed. There were studies which showed that smaller firms grow faster  or slower  than bigger ones. However, it is generally thought that the proportional rate of growth of a firm (conditional on survival) is decreasing in size, as far as small and medium firms are concerned, which share a large fraction of industrial sectors in number. However, the remaining larger firms constitute a small fraction in number, but occupy a large fraction of total sum of firms size. This is due to the effect of heavy tail, much heavier than expected from log-normal regime. See recent works 333 observed that log-normal distribution overestimates the upper-tail of size distribution based on Computat in U.S. As noted in the paper, the dataset is consisting of only publicly-traded firms. This can be a possible cause of their observation.  used much larger dataset in U.K. Though their plot showed a power-law regime over several orders of magnitude, they rejected the hypothesis of power-law due to the presence of super-giant firms. We consider that both of the these points deserve further investigation.. This is the Pareto-Zipf regime which we focus on in this paper.

On the other hand, the assumption (b) about temporal correlation between successive growth rates are not investigated with definite conclusions. , for example, showed that the distribution of growth rates shows a first-order positive autocorrelation: the growth process will result faster for firms which recorded a sharp growth in previous years. ( also furnished a test test for the validity of Gibrat’s Law that takes into account the “historical memory” of the growth process.)

Gibrat’s work also opened up a stream of theoretical models and ideas. Kalecki noted that Gibrat’s model leads to “unrealistic” feature, that is, the variance of the size distribution would increase indefinitely with time. He considered several models, one of which assumed that the expected rate of growth increased less than proportionately, leading to a log-normal distribution with constant variance.

Herbert Simon considered it more important that the firm size distribution has heavy tail in upper-region of size, which was better fitted by Yule distribution or asymptotically a Pareto-Zipf law. In order to explain such a distribution, based on his earlier work  for the explanation of Zipf’s law in word frequency, he assumed Gibrat’s law (in a much weaker form than ours) with a boundary condition for entry and exit of firms. In conformity with preceding work by Champernowne for personal income , Simon could show that the emergence of power-law behavior is quite robust irrespectively of modification of the stochastic process (see  for collection of related papers). Simon modeled the process of entry corresponding to new firms which compete with existing firms to catch market opportunities. This line of models was followed by  which relaxes the assumption of Gibrat’s law, and also by  which explained the Laplace distribution for growth rate. Simon also extended his model incorporating merger and acquisition process (see also  for recent work). These works attempted to take into account the direct and indirect interactions among firms, which was ignored in the assumption (c) above.

Our work is in affinity with Simon’s view in the points that the upper-tail of size distribution, Pareto-Zipf law, is focused rather than the log-normal regime, and that the origin of it is related to Gibrat’s law and boundary condition of entry-exit of firms. It is interesting to point out that Mansfield , following Simon’s model, empirically showed that the Gibrat’s law seemed to hold only above a certain minimum size of firms. (See also  for the influence of  onto later work.)

At the end of this brief survey, let us point out why recent advent of econophysics can have important impact on economics. The econophysics approach attempts to treat the whole industrial organization as a complex system, in which firms are interacting atoms, that exhibits universal scaling laws .

Concerning firm size, the Pareto-Zipf power-law distribution has a long history since the seminal work by Herbert Simon, but its study extending to the details of growth rate was only recently facilitated by modern datasets with good abundance and quality. In this line of research, resent findings (e.g. ) showed that power-law distribution gives a very good fit for different samples of firm size. In this paper we shall not only confirm this fact with different European countries and for different measures of size, but also uncover the underlying kinematics that relates Pareto-Zipf law to Gibrat law explicitly. Following the notion of self-organized criticality , the occurrence of a power-law reveals that a deep interaction among system’s subunits, reacting to idiosyncratic shocks, leads to a critical state in which no attractive points nor states emerge. Such interaction and critical states are so important notions with that of self-organized criticality. Under economic point of view, interaction means that it is not possible to define a representative agent because the dynamics of the system is originated just from the interaction among heterogeneous agents. Moreover, in consequence of critical state, equilibrium exists only as asymptote, along which the system moves from an unstable critical point to another. The authors believe that economics can enjoy these ideas coming from econophysics on heterogeneous interacting agents (see  for an example).

## 3 Dataset of European Firms

We use the dataset, Bureau van Dijk’s AMADEUS, which contains descriptive and balance data of about 260,000 firms of 45 European countries for the years 1992–2001. For every firm are reported a number of juridical, historical and descriptive data (as e.g. year of inclusion, participations, mergers and acquisitions, names of the board directors, news, etc.) and a series of data drawn from its balance and normalized. It reports the current values (for several currencies) of stocktaking, balance sheet (BS), profit and loss account (P/L) and ratios. The descriptive data are frequently updated while the numerical ones are taken from the last available balance. Since balance year does not always match conventional year, the number of firms included may vary during the year if one of the excluded firms in last recording satisfy one of the criteria described below. The amount and the completeness of available data differs from country to country. To be included in the data set firms must satisfy at least one of these three dimensional criteria:

* for U.K., France, Germany, Italy, Russian Federation and Ukraine

+ operating revenue equal to at least 15 million Euro
    + total assets equal to at least 30 million Euro
    + number of employees equal to at least 150
* for the other countries

+ operating revenue equal to at least 10 million Euro
    + total assets equal to at least 20 million Euro
    + number of employees equal to at least 100

As a proxy for firm size, we utilize one of the financial and fundamental variables; total-assets, sales and number of employees. We use number of employees as a complementary variable so as to check the validity and robustness of our results. Note that the dataset includes firms with smaller total-assets, simply because either the number of employees or the operating revenue (or both of them) exceeds the corresponding threshold. We thus focus on complete sets of those firms that have larger amount of total-assets than the threshold, and similarly those for number of employees. For sales, we assume that our dataset is nearly complete since a firm with a small amount of total-assets and a small number of employees is unlikely to make a large amount of sales. For our purposes, therefore, we discard all the data below each corresponding threshold for each measure of firm size. This procedure makes the number of data points much less. However, for a several developed countries, we have enough amount of data for the study of Gibrat’s law. In what follows, our results are shown for UK and France, although we obtained similar results for other developed countries. The threshold for total-assets in these two countries is 30 million euros, and that for number of employees is 150 persons, as described above. For sales, we used 15 million euros per year as a threshold. We will also show results for Italy and Spain in addition to U.K. and France only when examining the annual change of Pareto indices.

It should be remarked that other problems in treating these data takes origin from the omission, in the on-line dataset of AMADEUS, of the date of upgrade, so that it is often not clear when a firm changed its juridical status, or went bankrupted or inactive. For some countries the indication activity/inactivity is not shown at all, so that it was impossible, even indirectly, to individuate the year of exit. Therefore, our study should be taken as the analysis conditional on survival of firms.

## 4 Firm Growth

In this section, our results are shown for UK and France, and for total-assets, number of employees and sales. Each list of firms is exhaustive in the way we described in the preceding section.

### 4.1 Pareto-Zipf distribution

First we show that the distribution of firm size obeys a power-law in the range of our observation whatever we take as a variable for firm size. Fig. 1 depicts the cumulative distributions for total-assets in France (a), sales in France (b), and number of employees in UK (c). The number of data points is respectively (a) 8313, (b) 15776 and (c) 15055.

Pareto-Zipf law states that the cumulative distribution $`P_>(x)`$ for firm size $`x`$ follows eq. (1). The power-law fit for $`xx_0`$, where $`x_0`$ denotes the threshold mentioned above for each measure of firm size, gives the values of $`\mu `$; (a) 0.886$`\pm `$0.005, (b) 0.896$`\pm `$0.011, (c) 0.995$`\pm `$0.013 (standard error at 99% significance level). $`\mu `$ is close to unity. Note that the power-law fit is quite well nearly three orders of magnitude in size of firms.

Pareto index is surprising stable in its value. Fig. 2 is a panel for the annual change of Pareto indices for four countries, Italy, Spain, France and U.K. estimated from total-assets, number of employees and sales (except U.K.). Different measures of firm size give reasonably same behavior. It is observed that the value $`\mu `$ is quite stable being close to unity in all the countries.

### 4.2 Gibrat’s law

Let us denote a firm’s size by $`x`$ and its two values at two successive points in time (i.e., two consecutive years) by $`x_1`$ and $`x_2`$. Growth rate is given by $`Rx_2/x_1`$. We also express the rate in terms of its logarithm, $`r\mathrm{log}_{10}R`$. We examine the probability density for the growth rate $`P(r|x_1)`$ on the condition that the firm size $`x_1`$ in an initial year is fixed.

For the conditioning we divide the range of $`x_1`$ into logarithmically equal bins. For the total-assets in the dataset (Fig. 3 (a)), the bins are taken as $`x_13\times [10^{7+0.4(n1)},10^{7+0.4n}]`$ (euros) with $`n=1,\mathrm{},5`$. For the sales in (b), $`x_11.5\times [10^{7+0.4(n1)},10^{7+0.4n}]`$ (euros) with $`n=1,\mathrm{},5`$. For the number of employees in (c), $`x_11.5\times [10^{2+0.4(n1)},10^{2+0.4n}]`$ (persons) with $`n=1,\mathrm{},5`$. In all the cases, the range of conditioning covers two orders of magnitude in each variable. We calculated the probability density function for $`r`$ for each bin, and checked the statistical dependence on $`x_1`$ by graphical method.

Fig. 3 is the probability density function $`P(r|x_1)`$ for each case. It should be noted that due to the limit $`x_1>x_0`$ and $`x_2>x_0`$, the data for large negative growth are not available. In all the cases, it is obvious that the function $`P(r|x_1)`$ has little statistical dependence on $`x_1`$, since all the curves for different $`n`$ collapse on a single curve. This means that the growth rate is independent of firm size in the initial year. That is, Gibrat’s law holds.

### 4.3 Time-reversal symmetry

The validity of Gibrat’s law in the Pareto-Zipf regime appears to be in disagreement with recent literature on firm growth. In the next section, we will show that this is not actually the case by proving that Gibrat and Pareto-Zipf are equivalent under an assumption. The assumption is detailed balance, whose validity is checked here.

Let us denote the joint probability distribution function for the variable $`x_1`$ and $`x_2`$ by $`P_{12}(x_1,x_2)`$. The detailed balance, or what we call time-reversal symmetry, is the assumption that $`P_{12}(x_1,x_2)=P_{12}(x_2,x_1)`$. The joint probabilities for our datasets are depicted in Fig. 4 as scatter-plots of individual firms.

We used two different methods to check the validity of time-reversal symmetry. One is an indirect way to check a non-trivial relationship between the growth-rate in positive side ($`r>0`$) and that in negative ($`r<0`$). That is, as we shall prove in the next section, the probability density distribution in positive and negative growth rates must satisfy the relation given by eq. (27), if the property of time-reversal symmetry holds. We fitted the cumulative distribution only for positive growth rate by a non-linear function, converted to density function, and predicted the form of distribution for negative growth rate by eq. (27) so as to compare with the actual observation (see Appendix for details). In each plot of Fig. 3, a solid line in the $`r>0`$ side is such a fit, and a broken line in the $`r<0`$ side is our prediction. The agreement with the actual observation is quite satisfactory, thereby supporting the validity of time-reversal symmetry.

The other way we took is a direct statistical test for the symmetry in the two arguments of $`P_{12}(x_1,x_2)`$. This can be done by two-dimensional Kolmogorov-Smirnov (K-S) test, which is not widely known but was developed by astrophysicists . This statistical test is not strictly non-parametric (like the well-known one-dimensional K-S test), but has little dependence on parent distribution except through coefficient of correlation. We compare the scatter-plot sample for $`P_{12}(x_1,x_2)`$ with another sample for $`x_1`$ and $`x_2`$ interchanged by making the null hypothesis that these two samples are taken from a same parent distribution. We used the logarithms $`\xi _1=\mathrm{log}x_1`$ and $`\xi _2=\mathrm{log}x_2`$, and added constants to $`\xi _1`$ and $`\xi _2`$ so that the average growth rate is zero. This addition (or multiplication in $`x_1`$ and $`x_2`$) is simply subtracting the nominal effects due to inflation, etc. We applied two-dimensional K-S test to the resulting samples. The null hypothesis is not rejected in 95% significance level in all the cases we studied.

## 5 Pareto-Zipf’s law and Gibrat’s law under detailed balance

In the preceding section, we have shown that both of Pareto-Zipf and Gibrat’s laws hold for large firms. This suggests that these two laws are closely related with each other. We show in this section that in fact they are equivalent to each other under the condition of detailed balance.

Let $`x`$ be a firm’s size, and let its two values at two successive points in time (i.e., two consecutive years) be denoted by $`x_1`$ and $`x_2`$. We denote the joint probability distribution function (pdf) for the variable $`x_1`$ and $`x_2`$ by $`P_{12}(x_1,x_2)`$. The joint pdf of $`x_1`$ and the growth rate $`R=x_2/x_1`$ is denoted by $`P_{1R}(x_1,R)`$. Since $`P_{12}(x_1,x_2)dx_1dx_2=P_{1R}(x_1,R)dx_1dR`$ under the change of variables from $`(x_1,x_2)`$ to $`(x_1,R)`$, these two pdf’s are related to each other as follows:

$$P_{1R}(x_1,\frac{x_2}{x_1})=x_1P_{12}(x_1,x_2).$$
(3)
We define conditional probabilities:

$`P_{1R}(x_1,R)`$ $`=`$ $`P_1(x_1)Q(R|x_1)`$ (4)
$`=`$ $`P_R(R)S(x_1|R),`$ (5)
Both $`P_1(x_1)`$ and $`P_R(R)`$ are marginal:

$`P_1(x_1)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}P_{1R}(x_1,R)𝑑R\left(={\displaystyle _0^{\mathrm{}}}P_{12}(x_1,x_2)𝑑x_2\right),`$ (6)
$`P_R(R)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}P_{1R}(x_1,R)𝑑x_1,`$ (7)
since the following normalizability conditions are satisfied:

$`1`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}Q(R|x_1)𝑑R,`$ (8)
$`1`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}S(x_1|R)𝑑x_1.`$ (9)
Three phenomenological properties can be summarized as follows.

1. Detailed Balance (Time-reversal symmetry):
The joint pdf $`P_{12}(x_1,x_2)`$ is a symmetric function:

$$P_{12}(x_1,x_2)=P_{12}(x_2,x_1).$$
(10)
2. Pareto-Zipf’s law:
The pdf $`P_1(x)`$ obeys power-law for large $`x`$:

$$P_1(x)x^{\mu 1},$$
(11)
for $`x\mathrm{}`$ with $`\mu >0`$.
3. Gibrat’s law:
The conditional probability $`Q(R|x)`$ is independent of $`x`$:

$$Q(R|x)=Q(R).$$
(12)
We note here that this holds only for large $`x`$, because we confirmed it in actual data only in that region, and because otherwise it leads to an inconsistency, as we will see shortly. This relation was called Universality in . All the arguments below is restricted in this region.

Before starting our discussion of interrelation between these properties, let us first rewrite the detailed balance condition (A) in terms of $`P_{1R}(x_1,R)`$:

$`P_{1R}(x_1,R)`$ $`=`$ $`x_1P_{12}(x_1,x_2)`$ (13)
$`=`$ $`x_1P_{12}(x_2,x_1)`$
$`=`$ $`{\displaystyle \frac{x_1}{x_2}}x_2P_{12}(x_2,x_1)`$
$`=`$ $`R^1P_{1R}(x_2,R^1),`$
where eq. (10) was used in the second line, and eq. (3) was used in the first and the third line. The above relation may be rewritten as follows by the use of the conditional probability $`Q(R|x_1)`$ in eq. (5);

$`{\displaystyle \frac{Q(R^1|x_2)}{Q(R|x_1)}}=R{\displaystyle \frac{P_1(x_1)}{P_1(x_2)}}.`$ (14)
In passing, it should be noted that eq. (13) leads to the following:

$`P_R(R)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}P_{1R}(x_1,R)𝑑x_1`$ (15)
$`=`$ $`{\displaystyle _0^{\mathrm{}}}R^1P_{1R}(x_2,R^1)𝑑x_1`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}R^2P_{1R}(x_2,R^1)𝑑x_2`$
$`=`$ $`R^2P_R\left(R^1\right)`$
where eq. (13) was used in the second line, and the third line is merely change of integration variable. This relation between the marginal growth-rate pdf $`P_R(R)`$ for positive growth ($`R>1`$) and negative growth ($`R<1`$) leads to the following relation, as it should:

$$_1^{\mathrm{}}P_R(R)𝑑R=_0^1P_R(R)𝑑R.$$
(16)
### 5.1 (A)+(C)$``$(B)

Let us first prove that the properties (A) and (C) lead to (B). By substituting the Gibrat’s law eq. (12) in eq. (14), we find the following:

$$\frac{P_1(x_1)}{P_1(x_2)}=\frac{1}{R}\frac{Q(R^1)}{Q(R)}.$$
(17)
This relation can be satisfied only by a power-law function eq. (11).

\[Proof\]
Let us rewrite eq. (17) as follows:

$$P_1(x)=G(R)P_1(Rx),$$
(18)
where $`x`$ denotes $`x_1`$, and $`G(R)`$ denotes the right-hand side of eq. (17), i.e.

$$G(R)\frac{1}{R}\frac{Q(R^1)}{Q(R)}.$$
(19)
We expand this equation around $`R=1`$ by denoting $`R=1+ϵ`$ with $`ϵ1`$ as

$`P_1(x)`$ $`=`$ $`G(1+ϵ)P_1((1+ϵ)x)`$ (20)
$`=`$ $`(1+G^{}(1)ϵ+\mathrm{})(P_1(x)+P_1^{}(x)ϵx+\mathrm{})`$
$`=`$ $`P_1(x)+ϵ(G^{}(1)P_1(x)+xP_1^{}(x))+O(ϵ^2),`$
where we used the fact that $`G(1)=1`$. We also assumed that the derivatives $`G^{}(1)`$ and $`P_1^{}(x)`$ exists in the above, whose validity should be checked against the results. From the above, we find that the following should be satisfied

$$G^{}(1)P_1(x)+xP_1^{}(x)=0,$$
(21)
whose solution is given by

$$P_1(x)=Cx^{G^{}(1)}.$$
(22)
This is the desired result, Pareto-Zipf’s law, and is consistent with the assumption made earlier that $`P_1^{}(x)`$ exists. By substituting the result eq. (22) in eq. (19) and eq. (17), we find that

$$G(R)=R^{G^{}(1)},$$
(23)
which is consistent with the assumption that $`G^{}(1)`$ exists.

\[Q.E.D.\]

From eq. (19) we may calculate $`G^{}(1)`$ in terms of derivatives of $`Q(R)`$. It should, however, be noted that $`Q(R)`$ has a cusp at $`R=1`$ as is apparent in Fig. 3, and therefore $`Q^{}(R)`$ is expected not to be continuous at $`R=1`$. Bearing this in mind, we calculate $`G(1+ϵ)`$ for $`0<ϵ1`$ as follows:

$`G(1+ϵ)`$ $``$ $`{\displaystyle \frac{1}{1+ϵ}}{\displaystyle \frac{Q(1ϵ)}{Q(1+ϵ)}}`$ (24)
$``$ $`(1ϵ){\displaystyle \frac{Q(1)ϵQ^{}(1)}{Q(1)+ϵQ^+(1)}}`$
$``$ $`G(1)+ϵ\left(1{\displaystyle \frac{Q^+(1)+Q^{}(1)}{Q(1)}}\right),`$
where we denoted the right-derivative and left-derivative of $`Q(R)`$ at $`R=1`$ by the signs $`+`$ and $``$ in the superscript, respectively. From the above, we find that

$$G^{}(1)=1\frac{Q^+(1)+Q^{}(1)}{Q(1)},$$
(25)
From eq. (22) and eq. (25), we find that

$$\frac{Q^+(1)+Q^{}(1)}{Q(1)}=\mu 2.$$
(26)
From eqs.(19) and (23), we find the following relation:

$$Q(R)=R^{\mu 2}Q(R^1),$$
(27)
which should be in contrast to eq. (15). This is related to the point that we mentioned in eq. (12): If the Gibrat’s law eq. (12) holds for all $`x[0,\mathrm{}]`$, then $`P_R(R)=Q(R)`$ from eq. (7). If so, eq. (27) contradicts to eq. (15) since $`\mu >0`$. Besides, the Pareto-Zipf’s law we derived from Gibrat’s law is not normalizable if it holds for any $`x`$. Therefore, Gibrat’s law should hold only for large $`x`$.

The result eq. (27) shows that the function $`Q(R)`$ is continuous at $`R=1`$, as is easily seen by substituting $`R=1+ϵ`$ with $`ϵ>0`$ on both hand side and taking the limit $`ϵ+0`$. Also, by taking the derivative of the both hand side and taking the limit in a similar manner, we can reproduce eq. (26).

### 5.2 (A)+(B)$``$ ?

Let us next examine what we obtain if we had only Pareto-Zipf’s law instead of Gibrat’s law under the detailed balance.

In this case, substituting the Pareto-Zipf’s law eq. (11) into eq. (14) we find that

$$\frac{Q(R^1|Rx)}{Q(R|x)}=R^{\mu +2},$$
(28)
where we denote $`x_1`$ by $`x`$ and $`x_2`$ by $`Rx`$. We now define a function $`H(z,x)`$ as

$$Q(R|x)=x^{\mu +2}H(R^{1/2}x,x).$$
(29)
It should be noted that this does not constrain $`Q(R|x)`$ in any way: arbitrary function of the variable $`R`$ and $`x`$ can be written in the form of eq. (29). By substituting eq. (29) into eq. (28), we find that

$$H(R^{1/2}x,Rx)=H(R^{1/2}x,x),$$
(30)
which means that the function $`H(z,x)`$ has the following invariance property.

$$H(z,x)=H(z,z^2/x).$$
(31)
Other than this constraint and some trivial constraint such as continuity, there is no nontrivial constraint on $`H(z,x)`$ or $`Q(R|x)`$.

The results eqs. (29) and (31) is a generalization of the property eq. (27) we found earlier . In fact, the property eq. (27) follows from eq. (31) in the special case:

$$H(z,x)=Q((z/x)^2)x^{\mu 2},$$
(32)
for which eq. (29) becomes eq. (12), namely the statement of Gibrat’s law.

### 5.3 (B)+(C)$``$(A)?

Let us discuss the last question: Under Pareto’s and Gibrat’s laws, what can we say about the detailed balance? In order to answer this, we use eq. (11) and eq. (12) to write $`P_{1R}(x,R)`$ for large $`x`$ as follows:

$$P_{1R}(x,R)=Ax^{\mu 1}Q(R),$$
(33)
where $`A`$ is a proportionality constant. According to eq. (13), the detailed balance is satisfied if this is equal to

$$R^1P_{1R}(xR,R^1)=Ax^{\mu 1}R^{\mu 2}Q(R^1),$$
(34)
where we used eq. (33). Therefore, we find that the detailed balance condition is equivalent to eq. (27) in this case.

Summarizing this section, we have proved that under the condition of detailed balance (A), if the Pareto-Zipf law (B) holds in a region of firm size, then the Gibrat’s law (C) must hold in the region, and vice versa. The condition (A) means detailed-balance. On the other hand, if both of (B) and (C) hold, (A) follows provided that eq. (27) holds. eq. (27) is our prediction which gives a non-trivial relation between positive growth ($`R>1`$) and negative ($`R<1`$). This kinematic relation was empirically verified in Fig. 3. See also previous work  for the validity of this relation in personal income and firms tax-income in Japan.

## 6 Summary

The distribution of firm size is quite often dominated by power-law in the upper tail over several orders of magnitude. This regime of Pareto-Zipf law is different from log-normal distribution in the lower and sometimes wider regime of firm size. The upper tail is occupied by a small number of firms, but they dominate a large fraction of total sum of firm size.

By using exhaustive datasets of those large firms and with different measures of firm size in Europe, we show that the Pareto-Zipf law holds as in eq. (1) for firm size $`x`$ larger than observational threshold $`x_0`$, and that Gibrat’s law of proportionate effect holds as in eq. (2) for successive sizes $`x_1`$ and $`x_2`$ exceeding $`x_0`$, stating that the growth rate of each firm is independent of initial size. We also find that detailed balance holds which means that the frequency of transition from $`x_1`$ to $`x_2`$ is statistically the same as that for its reverse process. The Gibrat’s law, Pareto-Zipf’s law and detailed balance condition are related to each other. We prove various relationships among them. It follows as one of the consequences that there exists a relation between the positive and negative sides of the distribution of growth rate via the Pareto index. The relation is confirmed empirically in our dataset of European firms.

We thank Edmondo Di Giuseppe for helping the preparation of raw datasets before analysis. Y. F. and W. S. thank Katsunori Shimohara for encouragement. Supported in part by grants from the Telecommunications Advancement Organization in Japan and from Japan Association for Cultural Exchange.

## Appendix A Fitting distribution of growth rate

For the purpose of fitting probability density function of positive growth rate ($`R>1`$), we used cumulative distribution of positive growth rate, defined by

$$P_>^+(R)=\text{Prob}\{𝑹>R|R>1\}.$$

$`P_>^+(R)`$ can be estimated, as usual, by size versus rank plot restricted only for $`R>1`$ as follows. Let the number of all firms with $`R>1`$ be $`N_+`$, and sort their growth rates in descending order: $`R^{(1)}R^{(2)}\mathrm{}R^{(k)}\mathrm{}R^{(N_+)}`$. Then the estimate is given by

$$P_>^+(R)=\frac{k}{N_+}=C^1_{V^{(k)}}P_{1R}(x_1,R)𝑑x_1𝑑R,$$

where $`V^{(k)}=\{(x_1,R)|x_1x_0,RR^{(k)}\}`$ ($`x_0`$ is the observational threshold mentioned in section 4.1), and $`C`$ is the normalization:

$$C=_{x_0}^{\mathrm{}}𝑑x_1_1^{\mathrm{}}𝑑RP_{1R}(x_1,R).$$

Using the observational fact that eq. (12) holds in the region $`\{(x_1,R)|x_1x_0,R1\}`$, the above equation for $`P_>^+(R)`$ reads

$$P_>^+(R)=C^1_{x_0}^{\mathrm{}}𝑑x_1P_1(x_1)𝑑x_1_R^{\mathrm{}}𝑑R^{}Q(R^{})=Q_0{}_{}{}^{1}_{R}^{\mathrm{}}Q(R^{})𝑑R^{},$$
(35)
where the normalization factor is written by

$$Q_0=_1^{\mathrm{}}Q(R^{})𝑑R^{}.$$
(36)
By taking derivative of eq. (35) with respect to $`R`$, it follows that

$$Q(R)=Q_0\frac{d}{dR}P_>^+(R).$$
(37)
We empirically found that the rank-size plot can be well fitted by a non-linear function of the form:

$$\mathrm{log}_{10}P_>^+(R=10^r)=a(1e^{br})crF(r),$$
(38)
where $`a`$, $`b`$ and $`c`$ are parameters. An example is given in Fig. 5 for France total-assets (2001/2000). Cumulative probabilities $`P_>^+(r|x_1)`$ (the left-hand side of eq. (38)) conditioned on an initial year’s total-assets are shown for each of the same bins used in Fig. 3 (a), but restricted to the data with positive $`r`$. The non-linear fit done by eq. (38) is represented by a solid and bold line in the figure. Note also that the curves for different bins almost collapse because of the statistical independence of $`x_1`$.

Under the change of variable, $`r=\mathrm{log}_{10}R`$, the probability density for $`r`$ defined by $`q(r)`$ is related to that for $`R`$ by

$$\mathrm{log}_{10}q(r)=\mathrm{log}_{10}Q(R=10^r)+r+\mathrm{log}_{10}(\mathrm{ln}10)$$
(39)
Therefore it follows from eq. (37) and eq. (39) that

$$\mathrm{log}_{10}q(r)=F(r)+\mathrm{log}_{10}\left[\frac{dF(r)}{dr}\right]+\mathrm{log}_{10}Q_0+\mathrm{log}_{10}(\mathrm{ln}10).$$
(40)
In each plot of Fig. 3, the solid curve is given by eq. (40), where $`P(r|x_1)`$ denotes the probability density function $`q(r)`$ for $`r`$, conditioned on an initial year’s size $`x_1`$.

The relation eq. (27) for positive ($`R>1`$) and negative ($`R<1`$) growth rates can be written in terms of $`q(r)`$ as

$$\mathrm{log}_{10}q(r)=\mu r+\mathrm{log}_{10}q(r),$$
(41)
which is easily shown by eq. (39). In each plot of Fig. 3, the dotted curve for negative growth rate ($`r<0`$) is obtained from the solid curve for positive one ($`r>0`$) through the relation eq. (41).


text:
# Degeneration of linear systems through fat points on 𝐾3 surfaces

## 1. Introduction

In this paper we assume the ground field is the field of the complex numbers.

Let $`S`$ be a smooth projective generic $`K3`$ surface (i.e. $`\mathrm{Pic}(S)`$) and let $`H`$ be the generator of $`\mathrm{Pic}(S)`$.

Consider $`n`$ points in general position on $`S`$, to each one of them associate a natural number $`m_i`$ called the multiplicity of the point and let $`n_j`$ be the number of $`p_i`$ with multiplicity $`m_i`$.

For a linear system of curves in $`|dH|`$ with $`n_j`$ general base points of multiplicity $`m_j`$ for $`j=1\mathrm{}k`$, define its virtual dimension $`v`$ as $`dim|dH|n_im_i(m_i+1)/2`$ and its expected dimension by $`e=\mathrm{max}\{v,1\}`$. If the dimension of the linear system is $`l`$, then $`vel`$.

Observe that it is possible to have $`e<l`$, since the conditions imposed by the points may be dependent. In this case we say that the system is special.

Linear systems through general fat points on rational surfaces have been studied by many authors (see e.g. \[Gim89, Har85, CM98, CM01\]), but, as far as we know, on $`K3`$ surfaces, no results on the non-speciality of such systems are known.

In section 3 we develop a technique to degenerate a $`K3`$ surface and linear systems through fat points on $`K3`$ surfaces, this degeneration is similar to the degeneration of the plane used in \[BZ03\].

In the subsequent section, we use this degeneration to prove that, a homogeneous linear system of curves in $`|dH|`$ with $`n`$ general base points of multiplicity $`m`$ is non-special if all linear systems of curves in $`|dH|`$ with one multiple base point are non-special.

Finally, in section 5, we prove conjecture \[DL03, Conjecture 2.3 (i)\] for homogeneous linear systems of curves in $`|dH|`$ with $`4^u9^w`$ fat points of multiplicity $`m`$, if either $`v1`$ or $`v1`$ and $`u>0`$ or $`2d1mod3`$.

## 2. Preliminaries

Let $`S`$ be a smooth projective generic $`K3`$ surface and let $`H`$ be the generator of $`\mathrm{Pic}(S)`$, then $`H`$ is ample, $`H^2=2g22`$ and $`h^0(H)=g+1`$; moreover $`H`$ is very ample if $`g3`$ and if $`g=2`$ $`H`$ defines a double covering of $`^2`$ branched at an irreducible sextic (see \[May72, Proposition 3\]).

Consider $`Q_1,\mathrm{},Q_n`$ points in general position on $`S`$, for each one of these points fix a multiplicity $`m_1,\mathrm{},m_n`$. Define $`=^\gamma (d,m_1,\mathrm{},m_r)`$, with $`\gamma =H^2=2g2`$, as the linear system of curves in $`|dH|`$ with multiplicity $`m_i`$ at $`Q_i`$ for all $`i`$. Let $`v`$ denote $`\mathrm{vdim}()`$, then, using $`dim|dH|=\frac{\gamma d^2}{2}+1`$, we obtain that

$$v=\frac{\gamma d^2}{2}+1\underset{i=1}{\overset{n}{}}\frac{m_i(m_i+1)}{2}$$
(2.1)
Let $`S^{}`$ be the blowing-up of $`S`$ along the points $`Q_1,\mathrm{},Q_n`$, $`\pi :S^{}S`$ the projection map and $`E_i`$ the exceptional divisor corresponding to $`Q_i`$. The linear system $``$ then corresponds to the system $`\pi ^{}(|dH|)𝒪_S^{}(m_1E_1\mathrm{}m_nE_n)`$. By abuse of notation, we will denote this linear system on the blowing-up also by $`^\gamma (d,m_1,\mathrm{},m_r)`$.

Similary, by $`^\gamma (d,m_1^{n_1},\mathrm{},m_y^{n_y})`$, we denote the linear system of curves in $`|dH|`$ with $`n_i`$ points (in general position) of multiplicity $`m_i`$ (for all $`i=1,\mathrm{},n`$) as well as the corresponding linear system on the blowing-up of $`S`$ along those $`n_1+\mathrm{}+n_y`$ general points.

Let $`Z=_{i=1}^nm_iQ_i`$ be the $`0`$-dimensional scheme defined by the multiple points and consider the exact sequence of sheaves:

where $`_Z`$ is the ideal sheaf of $`Z`$. Taking cohomology we obtain

$$v=h^0(𝒪_S(dH)(Z))h^1(𝒪_S(dH)(Z))1,$$
(2.2)
because $`h^1(𝒪_S(dH))=0`$ (see e.g. \[May72\]).

Let $`e`$ denote the expected dimension of the linear system $``$, i.e. $`e=\mathrm{max}\{1,v\}`$ and let $`l`$ denote the dimension of $``$. Then obviously we have that $`vel`$.

If $`e<l`$ (resp. $`e=l`$), we say the linear system $``$ is special (resp. non-special). Note that equation (2.2) shows that a non-empty linear system $``$ is non-special if and only if $`h^1()=0`$.

## 3. Degeneration of linear systems on $`K3`$ surfaces

### 3.1. The degeneration of a $`K3`$ surface

Let $`S`$ be a $`K3`$ surface and $`\mathrm{\Delta }`$ a complex disc around the origin. Consider the product $`V=S\times \mathrm{\Delta }`$ and its two projections $`q_1:VS`$ and $`q_2:V\mathrm{\Delta }`$. Let $`V_t`$ denote $`S\times \{t\}`$. Consider $`b`$ general points in $`V_0`$ and blow up $`V`$ along those $`b`$ points. We then get a new threefold $`X`$ and the maps $`\pi :XV`$, $`p_1=q_1\pi :X\mathrm{\Delta }`$ and $`p_2=q_2\pi :XS`$; i.e. we obtain the following commutative diagram

Let $`X_t`$ be the fiber of $`p_1`$ over $`t\mathrm{\Delta }`$. If $`t0`$ then $`X_tV_t`$ is our $`K3`$ surface $`S`$. But, $`X_0`$ is the union of the proper transform $`\stackrel{~}{S}`$ of $`V_0`$ and the $`b`$ exceptional divisors $`𝐏_i`$. Obviously each $`𝐏_i`$ is isomorphic to $`^2`$ and $`\stackrel{~}{S}`$ is the blowing-up of our $`K3`$ surface $`S`$ at the $`b`$ general points with projection map $`𝔟:\stackrel{~}{S}S`$.

Every $`𝐏_i`$ intersects $`\stackrel{~}{S}`$ transversally along a curve $`R_i`$, which is a line in $`𝐏_i`$ and an exceptional divisor on $`\stackrel{~}{S}`$. When we want to indicate that we consider $`R_i`$ in $`\stackrel{~}{S}`$, resp. $`𝐏_i`$, we denote it by $`E_i`$, resp. $`L_i`$.

Note that the map $`p_1`$ gives a flat family of surfaces over $`\mathrm{\Delta }`$, so $`X_0`$ can be seen as a degeneration of $`S`$.

### 3.2. The degeneration of a linear system on a $`K3`$ surface $`S`$

Let $``$ be a line bundle on $`S`$, and, for $`k`$, define the line bundle $`𝒪_X(,k)`$ on $`X`$ by

$$𝒪_X(,k)=p_2^{}()𝒪_X(k\stackrel{~}{S}).$$

The restriction of $`𝒪_X(,k)`$ to $`X_t`$, for $`t0`$, is then isomorphic to $``$.

The restriction of $`𝒪_X(,k)`$ to $`X_0`$, which we denote by $`𝒳(,k)`$, is a flat limit of the line bundle $``$ on the general fiber $`X_t`$, so $`𝒳(,k)`$ can be seen as a degeneration of the line bundle $``$.

On any $`𝐏_i`$, we have that $`𝒳(,k)|_{𝐏_i}=𝒪_X(,k)|_{𝐏_i}𝒪_^2(k)`$, since $`\stackrel{~}{S}`$ intersects $`𝐏_i`$ along a line.

On $`\stackrel{~}{S}`$ we obtain that

$$𝒳(,k)|_{\stackrel{~}{S}}=𝒪_X(,k)|_{\stackrel{~}{S}}=p_2^{}()|_{\stackrel{~}{S}}𝒪_X(k\stackrel{~}{S})|_{\stackrel{~}{S}}.$$

But, since $`\stackrel{~}{S}X_t_{i=1}^l𝐏_i`$ as divisors on $`X`$ and $`𝐏_i`$ intersects $`\stackrel{~}{S}`$ along $`R_i`$, this means that

$$𝒳(,k)|_{\stackrel{~}{S}}𝔟^{}()𝒪_{\stackrel{~}{S}}(\underset{i=1}{\overset{b}{}}kE_i).$$

Let $`Q_1,\mathrm{},Q_n`$ be general points on $`S`$, consider a zero-dimensional subscheme $`Z=m_1Q_1+\mathrm{}+m_nQ_n`$ and let $``$ be the sheaf $`_Z`$, where $`_Z`$ denotes the ideal sheaf that defines $`Z`$. Choose positive integers $`a_1,\mathrm{},a_b`$ such that $`a_1+\mathrm{}a_bn`$. Now, for all $`i\{1,\mathrm{},b\}`$, consider $`a_i`$ general points on $`𝐏_i`$; and take $`n_{i=1}^ba_i`$ general points on $`\stackrel{~}{S}`$. Denote those $`n`$ points (on the $`𝐏_i`$ and $`\stackrel{~}{S}`$) by $`Q_i^{}`$ (any order will do), and let $`Z^{}`$ be the zero-dimensional subscheme of $`X_0`$ given by $`m_1Q_1^{}+\mathrm{}+m_nQ_n^{}`$.

Then we obtain that $`𝒳(,k)_Z^{}`$ is a degeneration of $``$.

### 3.3. Homogeneous linear systems on generic $`K3`$ surfaces

Let $`S`$ be a generic $`K3`$ surface and consider a homogeneous linear system $`=^\gamma (d,m^n)`$ on $`S`$.

Choose positive integers $`b`$ and $`a`$ such that $`abn`$, let $`X`$ and $`𝒳(,k)`$ be as constructed before. For all $`1ib`$ and $`1ja`$, let $`Q_{i,j}^{}`$ be a general point on $`𝐏_i`$ and for $`1jnab`$ let $`Q_i^{}`$ be a general point on $`\stackrel{~}{S}`$. Now consider the zero-dimensional subscheme

$$Z^{}=\underset{i=1}{\overset{nab}{}}mQ_i^{}+\underset{\begin{array}{c}i=1,\mathrm{},b\\ j=1,\mathrm{},a\end{array}}{}mQ_{i,j}^{},$$

On $`X_0`$ define $`_0:=𝒳(,k)_Z^{}`$, $`_i:=_0|_{𝐏_i}`$ and $`_{\stackrel{~}{S}}:=_0|_{\stackrel{~}{S}}`$. Then we see that a divisor in the linear system $`_0`$ on $`X_0`$ consists of a divisor $`D_{\stackrel{~}{S}}_{\stackrel{~}{S}}`$ and divisors $`D_i_i`$ such that $`D_{\stackrel{~}{S}}|_{R_i}=D_i|_{R_i}`$ for all $`i`$.

###### Notation 3.1.

1.
2. $`l=dim`$, $`l_0=dim_0`$, $`l_𝐏=dim_i`$, $`l_{\stackrel{~}{S}}=dim_{\stackrel{~}{S}}`$
3. $`_{\stackrel{~}{S}}=_{\stackrel{~}{S}}|_{_{i=1}^bE_i}`$, $`r_{\stackrel{~}{S}}=dim_{\stackrel{~}{S}}`$
4. $`_i=_𝐏|_{L_i}`$, $`r_𝐏=dim_i`$
5. $`\widehat{}_{\stackrel{~}{S}}=_{\stackrel{~}{S}}𝒪_{\stackrel{~}{S}}(_{i=1}^bE_i)`$, $`\widehat{l}_{\stackrel{~}{S}}=dim\widehat{}_{\stackrel{~}{S}}`$
6. $`\widehat{}_i=_i𝒪_{𝐏_i}(L_i)`$, $`\widehat{l}_𝐏=dim\widehat{}_i`$

Obviously we have the following

$`l_{\stackrel{~}{S}}`$ $`=r_{\stackrel{~}{S}}+\widehat{l}_{\stackrel{~}{S}}+1`$ (3.1)
$`l_𝐏`$ $`=r_𝐏+\widehat{l}_𝐏+1.`$ (3.2)
Let $`_{L_i}`$ denote the linear system on $`_{i=1}^bR_i`$ for which $`_{L_i}|_{R_i}=_i`$ and denote $`r_{L_i}=dim_{L_i}`$. Then we have the following equality

$$l_0=dim(_{\stackrel{~}{S}}_{L_i})+b(\widehat{l}_𝐏+1)+\widehat{l}_{\stackrel{~}{S}}+1$$
(3.3)
But, proceeding as in \[BZ03, § 2\], we obtain the transversality of $`_{\stackrel{~}{S}}`$ and $`_{L_i}`$; i.e.

$`dim(_{\stackrel{~}{S}}_{L_i})`$ $`=\mathrm{max}\{1,r_{\stackrel{~}{S}}+r_{L_i}bk\}`$
$`=\mathrm{max}\{1,r_{\stackrel{~}{S}}+br_𝐏bk\}.`$ (3.4)
For the virtual dimensions of the systems we introduce the following

###### Notation 3.2.

1.
2. $`v=\mathrm{vdim}`$
3. $`v_𝐏=\mathrm{vdim}_i`$, $`v_{\stackrel{~}{S}}=\mathrm{vdim}_{\stackrel{~}{S}}`$
4. $`\widehat{v}_{\stackrel{~}{S}}=\mathrm{vdim}\widehat{}_{\stackrel{~}{S}}`$, $`\widehat{v}_𝐏=\mathrm{vdim}\widehat{}_i`$

Since $`=^\gamma (d,m^n)`$, $`_{\stackrel{~}{S}}^\gamma (d,k^b,m^{nab})`$, $`\widehat{}_{\stackrel{~}{S}}^\gamma (d,(k+1)^b,m^{nab})`$, $`_i(k,m^a)`$ and $`\widehat{}_i(k1,m^a)`$; a simple calculation shows that

$`v`$ $`=v_{\stackrel{~}{S}}+b\widehat{v}_𝐏+b=v_{\stackrel{~}{S}}+b(v_𝐏k)`$
$`=\widehat{v}_{\stackrel{~}{S}}+bv_𝐏+b=\widehat{v}_{\stackrel{~}{S}}+b(\widehat{v}_𝐏+k+2)`$ (3.5)
### 3.4. Remark

Let $`S`$ be a $`K3`$ surface, $``$ a line bundle on $`S`$ and let $`X`$ and $`𝒳(,k)`$ be as before.

Now do the construction of § 3.1 and § 3.2 using $`X_0`$ in stead of $`S`$ and $`𝒳(,k)`$ (or $`𝒳(,k)_Z^{}`$) in stead of $``$ (or $`_Z`$); i.e. consider $`W=X_0\times \mathrm{\Delta }`$, blow up $`W`$ along $`b^{}`$ general points on $`\stackrel{~}{S}\times \{0\}W_0=X_0\times \{0\}`$ and obtain the following commutative diagram

So we obtain a degeneration of $`X_0`$, and a degeneration of $`𝒳(,k)`$ (or $`𝒳(,k)_Z^{}`$).

We call this a double degeneration of $`S`$ and $``$ (or $`_Z`$), and, continuing in the same way, we can obtain an $`\eta `$-uple degeneration of $`S`$ and $``$ (or $`_Z`$), for any $`\eta 2`$.

## 4. Applying the degeneration to homogeneous linear systems on generic $`K3`$ surfaces

From now on we assume that we work on a generic $`K3`$ surface $`S`$, with $`H^2=\gamma `$, where $`H`$ is the generator of $`\mathrm{Pic}S`$.

The main result of this section is the following

###### Theorem 4.1.

If $`^\gamma (d,\mu )`$ is non-special for all $`\mu `$, then $`=^\gamma (d,m^n)`$ with $`n=4^u9^w`$ is non-special for all positive integers $`m,u`$ and $`w`$.

To make the proof of this theorem more transparent we first fix some notation and state a few auxiliary results.

Let $`c\{4,9\}`$ such that $`c|n`$ and consider the degeneration $`(X_0,_0)`$ of $`(S,)`$ obtained by using the construction explained in § 3 with $`b=n/c`$ and $`k`$. Note that, since $`_0`$ is a degeneration of $``$, $`l_0lev`$.

On $`^2`$, let $`(\delta ,\mu ^\nu )`$ denote the linear system of plane curves of degree $`\delta `$ having multiplicity $`\mu `$ at $`\nu `$ points (in general position).

The following is a well know result, and can easily be checked using for instance the results of \[Har85\].

###### Lemma 4.2.

On $`^2`$ the linear system $`(\delta ,\mu ^c)`$ with $`c\{4,9\}`$ is non-special for all $`\delta `$ and $`\mu `$. $`\mathrm{}`$

###### Remark 4.3.

Lemma 4.2 implies in particular that $`_i`$ and $`\widehat{}_i`$ are non-special linear systems on $`𝐏_i`$.

###### Claim 4.4.

If $`v1`$, $`v_{\stackrel{~}{S}}1`$, $`v_𝐏1`$ and $`_{\stackrel{~}{S}}`$ and $`\widehat{}_{\stackrel{~}{S}}`$ are non-special systems; then $`dim(_{\stackrel{~}{S}}_{L_i})=r_{\stackrel{~}{S}}+br_𝐏bk`$.

###### Proof.

Note that the conditions of the theorem immediately imply that $`l_{\stackrel{~}{S}}=v_{\stackrel{~}{S}}`$ and $`l_𝐏=v_𝐏`$.

Since $`dim(_{\stackrel{~}{S}}_{L_i})=\mathrm{max}\{1,r_{\stackrel{~}{S}}+br_𝐏bk\}`$, it suffices to show that $`r_{\stackrel{~}{S}}+br_𝐏bk1`$.

Because of (3.1) and (3.2), we have that

$$r_{\stackrel{~}{S}}+br_𝐏bk=l_{\stackrel{~}{S}}\widehat{l}_{\stackrel{~}{S}}1+b(l_𝐏\widehat{l}_𝐏1)bk.$$

If $`\widehat{v}_{\stackrel{~}{S}}1`$ and $`\widehat{v}_𝐏1`$, then $`\widehat{l}_{\stackrel{~}{S}}=\widehat{v}_{\stackrel{~}{S}}`$ and $`\widehat{l}_𝐏=\widehat{v}_𝐏`$. So

$`r_{\stackrel{~}{S}}+br_𝐏bk`$ $`=v_{\stackrel{~}{S}}\widehat{v}_{\stackrel{~}{S}}1+b(v_𝐏\widehat{v}_𝐏1)bk`$
$`=b(k+1)11.`$
If $`\widehat{v}_{\stackrel{~}{S}}2`$ and $`\widehat{v}_𝐏1`$, then $`\widehat{l}_{\stackrel{~}{S}}=1`$ and $`\widehat{l}_𝐏=\widehat{v}_𝐏`$. So

$`r_{\stackrel{~}{S}}+br_𝐏bk`$ $`=v_{\stackrel{~}{S}}+b(v_𝐏\widehat{v}_𝐏1)bk`$
$`=v_{\stackrel{~}{S}}1.`$
If $`\widehat{v}_{\stackrel{~}{S}}1`$ and $`\widehat{v}_𝐏2`$, then $`\widehat{l}_{\stackrel{~}{S}}=\widehat{v}_{\stackrel{~}{S}}`$ and $`\widehat{l}_𝐏=1`$. So

$`r_{\stackrel{~}{S}}+br_𝐏bk`$ $`=v_{\stackrel{~}{S}}\widehat{v}_{\stackrel{~}{S}}1+bv_𝐏bk`$
$`=b(1+v_𝐏)11.`$
If $`\widehat{v}_{\stackrel{~}{S}}2`$ and $`\widehat{v}_𝐏2`$, then $`\widehat{l}_{\stackrel{~}{S}}=1`$ and $`\widehat{l}_𝐏=1`$. So

$`r_{\stackrel{~}{S}}+br_𝐏bk`$ $`=v_{\stackrel{~}{S}}+b(v_𝐏k)`$
$`\stackrel{\text{(}\text{3.5}\text{)}}{=}v1.`$
###### Lemma 4.5.

If $`v1`$, $`v_{\stackrel{~}{S}}1`$, $`v_𝐏1`$ and $`_{\stackrel{~}{S}}`$ and $`\widehat{}_{\stackrel{~}{S}}`$ are non-special systems; then $``$ is non-special.

###### Proof.

Because $`_0`$ is a degeneration of $``$, we know that $`vll_0`$; so it suffices to prove that $`l_0=v`$.

Using (3.3), (3.4) and claim 4.4 we obtain that

$$l_0=r_{\stackrel{~}{S}}+br_𝐏bk+b(\widehat{l}_𝐏+1)+\widehat{l}_{\stackrel{~}{S}}+1.$$

So, using (3.1) and (3.2), we see that

$$l_0=l_{\stackrel{~}{S}}+b(l_𝐏k)=v_{\stackrel{~}{S}}+b(v_𝐏k)\stackrel{\text{(}\text{3.5}\text{)}}{=}v.$$

###### Lemma 4.6.

Let $`=^\gamma (d,m^n)`$ with $`n=4^u9^w`$, $`d,m,u,w`$, $`u+w>0`$ and $`v1`$. Take $`c`$, $`b`$, $`X_0`$ and $`_0`$ as before. Then $`k`$ such that $`v_{\stackrel{~}{S}}1`$ and $`v_𝐏1`$.

###### Proof.

Because $`v_{\stackrel{~}{S}}=\frac{\gamma }{2}d^2+1b\frac{k(k+1)}{2}`$ and $`v_𝐏=\frac{k(k+3)}{2}c\frac{m(m+1)}{2}`$, the inequalities $`v_{\stackrel{~}{S}}1`$ and $`v_𝐏1`$ are equivalent to

$`k^2+k\alpha 0,`$ $`\text{with }\alpha ={\displaystyle \frac{1}{b}}(\gamma d^2+4)\text{and}`$ (4.1)
$`k^2+3k\beta 0,`$ $`\text{with }\beta =cm(m+1)2.`$ (4.2)
But this is the same as

$`k`$ $`[{\displaystyle \frac{1\sqrt{1+4\alpha }}{2}},{\displaystyle \frac{1+\sqrt{1+4\alpha }}{2}}]\text{and}`$ (4.3)
$`k`$ $`]\mathrm{},{\displaystyle \frac{3\sqrt{9+4\beta }}{2}}][{\displaystyle \frac{3+\sqrt{9+4\beta }}{2}},\mathrm{}[.`$ (4.4)
So proving the statement is equivalent to proving that there exists a positive integer $`k`$ such that both (4.3) and (4.4) are satisfied, i.e. it is enough to show that

$$\frac{3+\sqrt{9+4\beta }}{2}+1\frac{1+\sqrt{1+4\alpha }}{2}.$$

A simple calculation shows that the previous inequality is equivalent to $`\alpha \beta +2`$, which is in turn equivalent to $`v1`$. ∎

###### Lemma 4.7.

If $`v1`$, $`\widehat{v}_{\stackrel{~}{S}}1`$, $`\widehat{v}_𝐏1`$ and $`_{\stackrel{~}{S}}`$ and $`\widehat{}_{\stackrel{~}{S}}`$ are non-special systems; then $``$ is non-special and thus empty.

###### Proof.

Since $`v_{\stackrel{~}{S}}=\widehat{v}_{\stackrel{~}{S}}+b(k+1)`$ and $`v_𝐏=\widehat{v}_𝐏+k+1`$, we see that $`v_{\stackrel{~}{S}}l(k+1)1`$ and $`v_𝐏k`$.

If $`v_{\stackrel{~}{S}}1`$ and $`v_𝐏1`$ then $`r_{\stackrel{~}{S}}=r_𝐏=1`$. If $`v_{\stackrel{~}{S}}1`$ and $`v_𝐏1`$, then $`r_{\stackrel{~}{S}}=1`$ and $`r_𝐏=v_𝐏k`$. If $`v_{\stackrel{~}{S}}1`$ and $`v_𝐏1`$, then $`r_{\stackrel{~}{S}}=v_{\stackrel{~}{S}}b(k+1)1`$ and $`r_𝐏=1`$. If $`v_{\stackrel{~}{S}}1`$ and $`v_𝐏1`$, then $`r_{\stackrel{~}{S}}=v_{\stackrel{~}{S}}`$ and $`r_𝐏=v_𝐏`$. So in any case we obtain that $`r_{\stackrel{~}{S}}+br_𝐏bk1`$, i.e. $`dim(_{\stackrel{~}{S}}_{L_i})=1`$.

Using (3.3) and $`\widehat{l}_{\stackrel{~}{S}}=\widehat{l}_𝐏=1`$, we see that $`l_0=\widehat{l}_{\stackrel{~}{S}}+k(\widehat{l}_𝐏+1)=1`$. ∎

###### Lemma 4.8.

Let $`=^\gamma (d,m^n)`$ with $`n=4^u9^w`$, $`d,m,u,w`$, $`u+w>0`$ and $`v1`$. Take $`c`$, $`b`$, $`X_0`$ and $`_0`$ as before. Then $`k`$ such that $`\widehat{v}_{\stackrel{~}{S}}1`$ and $`\widehat{v}_𝐏1`$.

###### Proof.

Because $`\widehat{v}_{\stackrel{~}{S}}=\frac{\gamma }{2}d^2+1b\frac{(k+1)(k+2)}{2}`$ and $`\widehat{v}_𝐏=\frac{(k1)(k+2)}{2}c\frac{m(m+1)}{2}`$, the inequalities $`\widehat{v}_{\stackrel{~}{S}}1`$ and $`\widehat{v}_𝐏1`$ are equivalent to

$`k^2+3k\alpha 0,`$ $`\text{with }\alpha ={\displaystyle \frac{1}{b}}(\gamma d^2+4)2\text{and}`$
$`k^2+k\beta 0,`$ $`\text{with }\beta =cm(m+1).`$
Proceeding as in the proof of lemma 4.6, we obtain that it is sufficient to prove that

$$\frac{1+\sqrt{1+4\beta }}{2}\frac{3+\sqrt{9+4\alpha }}{2}+1,$$

which is equivalent to $`\beta \alpha +2`$. And this last inequality is equivalent to $`v1`$. ∎

###### Proof of theorem 4.1.

Using induction, the result follows immediately from lemma’s 4.5 and 4.6 if $`v1`$, and from lemma’s 4.7 and 4.8 if $`v1`$. ∎

### 4.1. Using a higher order degeneration

Using the $`\eta `$-uple degeneration explained in Remark3.4, we can, proceeding as above, prove that the non-speciality of $`^\gamma (d,\mu _1,\mathrm{},\mu _\eta )`$ for all $`\mu _1,\mathrm{},\mu _\eta `$ implies the non-speciality of $`=^\gamma (d,m_1^{n_1},\mathrm{},m_\eta ^{n_\eta })`$ with $`n_i=4^{u_i}9^{w_i}`$ for all $`m_1,\mathrm{},m_\eta ,u_1,\mathrm{},u_\eta ,w_1,\mathrm{},w_\eta `$.

As we will see in section 5.3, proving the non-speciality of systems $`^\gamma (d,\mu )`$ is already rather complex. In fact, in section 5.3 we only prove when such a system is non-special for general generic $`K3`$ surfaces with $`\gamma =4`$.

### 4.2. Remark

As proved in \[CM01, Theorem 6.1\], the Segre conjecture on planar linear systems through fat points implies that the only special homogeneous linear systems $`(k,m^n)`$ are the ones with $`n\{2,3,5,6,7,8\}`$. So if the Segre conjecture is true, we can do the degeneration using, not only $`a=4`$ or $`9`$, but $`a\{4\}_9`$. In this way we can then prove theorem 4.1, for any $`n`$ which can be written as a product of powers of numbers in $`\{4\}_9`$, and thus in particular for any $`n9`$.

## 5. The non-speciality of linear systems $`^\gamma (d,m^n)`$ with $`n=4^u9^w`$ on general generic $`K3`$ surfaces with $`\gamma =4`$

The main result of this section is the following

###### Theorem 5.1.

Let $`S`$ be a general generic $`K3`$ surface with $`\gamma =4`$ and let $`=^4(d,m^n)`$, with $`n=4^u9^w`$, $`u,w_0`$ and $`d_{>0}`$.
(1) If $`v()1`$ then the linear system $``$ is non-special.
(2) If $`v()1`$ and either $`u>0`$ or $`2d1mod3`$ then $``$ is non-special unless $`n=1`$, $`m=2d`$ and $`d2`$. In the latter case $`dim(^4(d,2d))=0>\mathrm{edim}(^4(d,2d))=1`$ and $`^4(d,2d)=dC`$, with $`C`$ the unique element of $`^4(1,2)`$.

###### Remark 5.2.

Note that the previous theorem implies that \[DL03, Conjecture 2.3 (i)\] is true for homogeneous linear systems with $`n=4^u9^w`$ fat points on $`S`$ in the following cases
(i) $`v1`$, or
(ii) $`v1`$ and $`u>0`$, or
(iii) $`v1`$ and $`2d1mod3`$.

We will prove theorem 5.1 by using the degeneration as introduced in § 4 and the following

###### Proposition 5.3.

Let $`S`$ be a general generic $`K3`$ surface with $`\gamma =4`$. Then $`=^4(d,\mu )`$, with $`d,\mu _{>0}`$, is non-special unless $`\mu =2d`$ and $`d2`$. In the latter case $`dim(^4(d,2d))=0>\mathrm{vdim}(^4(d,2d))=1`$ and $`^4(d,2d)=dC`$, with $`C`$ the unique element of $`^4(1,2)`$.

###### Proof.

The map $`\varphi :S^3`$ corresponding to $`|H|`$ (with $`H`$ the generator of $`\mathrm{Pic}S`$) is an embedding, so we may look at $`S`$ as a quartic surface in $`^3`$. The unique element $`C`$ of $`^4(1,2)`$ is then the divisor on $`S^3`$ defined by the tangent plane to $`S`$ at this general point $`P`$. So $`C`$ is an irreducible plane curve of degree 4 having a node at $`P`$ (and no other singularities).

Blow up $`S`$ along the point $`P`$, and, by abuse of notation, let $`^4(d,\mu )`$ also denote the complete linear system on this blowing-up corresponding to $`^4(d,\mu )`$ on $`S`$. Let $`\stackrel{~}{C}`$ denote the strict transform of $`C`$ on this blowing-up, then $`g(\stackrel{~}{C})=2`$.

Assume $`\mu =2d`$ ($`d2`$) and consider the following sequence

Because $`\mathrm{deg}(^4(d,2d)𝒪_{\stackrel{~}{C}})=0=2g(\stackrel{~}{C})4`$, we know that $`h^0(^4(d,2d)𝒪_{\stackrel{~}{C}})=0`$ unless $`^4(d,2d)𝒪_{\stackrel{~}{C}}=𝒪_{\stackrel{~}{C}}`$. But $`^4(d,2d)𝒪_{\stackrel{~}{C}}=|dg_4^22d(P_1+P_2)|`$, where $`P_1`$ and $`P_2`$ are the intersection points of $`\stackrel{~}{C}`$ with the exceptional curve. So this would mean that $`|dg_4^22d(P_1+P_2)|=𝒪_{\stackrel{~}{C}}`$. But, since $`P`$ is a general point on the general quaritc $`S`$, we know that $`g_4^2=|K_{\stackrel{~}{C}}+P_1+P_2|`$; so we would obtain that $`|dK_{\stackrel{~}{C}}d(P_1+P_2)|=𝒪_{\stackrel{~}{C}}`$, or thus that $`|dK_{\stackrel{~}{C}}|=|d(P_1+P_2)|`$, which is not true (because $`P`$ is a general element of $`S`$). So $`h^0(^4(d,2d))=h^0(^4(d1,2(d1)))`$ for all $`d2`$. Using this a number of times we thus obtain that $`h^0(^4(d,2d))=h^0(^4(1,2))=1`$, so $`dC`$ is the only divisor in $`^4(d,2d)`$.

It now follows immediately that $`dim=1`$ if $`\mu 2d+1`$, since the only divisor in $`^4(d,2d)`$ has multiplicity exactly $`2d`$ in $`P`$.

Now consider the case $`\mu 2d1`$. If we can prove that $``$ is non-special for $`\mu =2d1`$, then the non-speciality follows for all $`\mu 2d1`$ (since $`\mathrm{vdim}^4(d,2d1)=d+1`$). For $`d=1`$, there is nothing to prove, since $`P`$ is a general point of $`S`$. So we may assume that $`d2`$ and that the $`^4(d^{},2d^{}1)`$ is non-special for $`d^{}d1`$. Now consider the following sequence

Because $`\mathrm{deg}(^4(d,2d1)𝒪_{\stackrel{~}{C}})=2=2g(\stackrel{~}{C})2`$, we know that $`h^1(^4(d,2d1)𝒪_{\stackrel{~}{C}})=0`$ unless $`^4(d,2d1)𝒪_{\stackrel{~}{C}}=K_{\stackrel{~}{C}}`$ (the canonical class on $`\stackrel{~}{C}`$). But $`^4(d,2d)𝒪_{\stackrel{~}{C}}=|dg_4^2(2d1)(P_1+P_2)|`$ and $`g_4^2=|K_{\stackrel{~}{C}}+P_1+P_2|`$ (where $`P_1`$ and $`P_2`$ are, as before, the intersection points of $`\stackrel{~}{C}`$ with the exceptional curve). Then this would mean that $`|(d1)K_{\stackrel{~}{C}}|=|(d1)(P_1+P_2)|`$, which is not true (because $`P`$ is a general element of $`S`$). So $`h^1(^4(d,2d1)𝒪_{\stackrel{~}{C}})=0`$ and, by hypotheses, $`h^1(^4(d1,2d3))=0`$, thus also $`h^1(^4(d,2d1))=0`$. ∎

###### Proof of theorem 5.1.

(1) In case $`n=1`$ we are done because of proposition 5.3. So we assume $`n>1`$, or equivalently $`u+w1`$. As in the proof of theorem 4.1, we can use lemma’s 4.5 and 4.6, but we have to end the induction with a class $`_{\stackrel{~}{S}}=^4(d,\mu )`$ such that $`\mu 2d`$ and $`\mu 2d1`$ (in order to apply lemma 4.5 we need to have $`_{\stackrel{~}{S}}`$ and $`\widehat{}_{\stackrel{~}{S}}`$ non-special). So, it is enough to prove that, in case $`u+w=1`$, we can take $`k`$ such that $`v_{\stackrel{~}{S}}1`$, $`v_𝐏1`$ and $`k\{2d,2d1\}`$. Since now $`v_{\stackrel{~}{S}}1`$, is equivalent to $`k2d1`$, it is enough to prove that $`v_𝐏1`$ for $`k=2d2`$ (which then implies that you can always take $`k=2d2`$ in the last step). Now, $`v_𝐏=\frac{k(k+3)}{2}c\frac{m(m+1)}{2}`$, so, for $`k=2d2`$, we obtain $`v_𝐏=2d^2d1c\frac{m(m+1)}{2}`$. Since $`u+w=1`$, we have that $`c=n`$, so $`c\frac{m(m+1)}{2}=2d^2+1v`$, and we see that $`v_𝐏1`$ for $`k=2d2`$ if and only if $`vd+1`$. Now, we only know that $`v1`$, and $`c=n\{4,9\}`$, but a simple calculation shows that in case $`n=4`$:

$$v1md1v1+2d.$$

In case $`n=9`$, we may assume that $`d4`$, since for $`d3`$, we obtain that $`v1`$ is only possible if $`m1`$. And, a simple calculation then shows that

$$v1m\frac{2d2}{3}v2+d$$

(note that, if $`m=\frac{2d1}{3}`$, we obtain $`v=2d<1`$, since $`d4`$).
(2) Again, if $`n=1`$ we are done because of proposition 5.3. So we assume $`n>1`$, or equivalently $`u+w1`$. As in the proof of theorem 4.1, we can use lemma’s 4.7 and 4.8, up to the last step. For this last step, we have $`c=n`$ and $`b=1`$, and we will take $`k=2d`$. In this case, we can no longer use lemma 4.7 since $`_{\stackrel{~}{S}}`$ is special. If we can prove however that $`dim(_{\stackrel{~}{S}}_{L_i})=\widehat{l}_{\stackrel{~}{S}}=\widehat{l}_𝐏=1`$, then we still obtain $`l_0=1`$, and thus also $`l=1`$. Since $`dim(_{\stackrel{~}{S}}_{L_i})=\mathrm{max}\{1,r_{\stackrel{~}{S}}+(r_𝐏2d)\}`$, $`r_{\stackrel{~}{S}}=0`$ and $`r_𝐏=l_𝐏`$ it is enough to prove that, for $`k=2d`$, $`v_𝐏2d1`$ and $`\widehat{v}_𝐏1`$. These inequalities are equivalent to respectively $`vd`$ and $`v1d`$, so it is sufficient to have $`vd`$.

In case $`u>0`$, we can make sure that in the last step $`c=4`$, and we obtain

$$v1mdv12d.$$

In case $`u=0`$, we have $`c=9`$ in the last step, and using the fact that $`2d1mod3`$, a simple calculation shows that

$$v1m\frac{2d}{3}v13d.$$

## Acknowledgement

Both authors would like to thank Marc Coppens for helpfull discussions on the subject.