The softening of wet lipid bilayer membranes during their gel-to-fluid first-order phase transition is studied by computer simulation of a family of two-dimensional microscopic interaction models. The models include a variable number, q, of lipid chain conformational states, where 2 less than or equal to q less than or equal to 10. Results are presented as functions of q and temperature for a number of bulk properties, such as internal energy, specific heat, and lateral compressibility. A quantitative account is given of the statistics of the lipid clusters which are found to form in the neighborhood of the transition. The occurrence of these clusters is related to the softening and the strong thermal density fluctuations which dominate the specific heat and the lateral compressibility for the high-q models. The cluster distributions and the fluctuations behave in a manner reminiscent of critical phenomena and percolation. The findings of long-lived metastable states and extremely slow relaxational behavior in the transition region are shown to be caused by the presence of intermediate lipid chain conformational states which kinetically stabilize the cluster distribution and the effective phase coexistence. This has as its macroscopic consequence that the first-order transition appears as a "continuous" transition, as invariably observed in all experiments on uncharged lecithin bilayer membranes. The results also suggest an explanation of the non-horizontal isotherms of lipid monolayers. Possible implications of lipid bilayer softening and enhanced passive permeability for the functioning of biological membranes are discussed.