Experimental and theoretical studies of the fluorescence intensity decays in biomolecular systems showed that under constraints of typical experiment fluorescence lifetime distribution is given by gamma function, which led to a power-like decay function. The latter well fits complex (heterogeneous) as well as simple mono-exponential decays. Fluorescence decay kinetics is described by mean value of lifetime distribution characterizing the average rate of the excited-state decay, and one new parameter of heterogeneity (1 < q < 3/2) describing the relative variance of distribution, and objectively reflecting physical heterogeneity of the system. In the classical limit, when q → 1, the gamma distribution becomes the Dirac delta function, and decay function converges from power-like form to the single-exponential form. Numerous examples illustrate successful applications of this model to rational analysis of complex fluorescence decays of biomacromolecules, e.g., complexes of E. coli purine nucleoside phosphorylase (PNP-I) with formycin A (FA, inhibitor) and orthophosphate (Pi, cosubstrate), which led to identification of the mechanism of deactivation of excited state and the N(2)H tautomeric form of FA selectively bound by PNP-I. The latter is of great importance for the studies of the mechanism of protein (enzyme) action as well as for more rational drug design.