function [fitness_gradient_f] = fitness_gradient_fish(F,a0F,a2ZM,GGE,mf,a2FM,mz,Z,h0F,W,a3FM,a0W,mw,a3WM,h0W,d0,dexp,h3WM,h3FM,h2FM,h2ZM) fitness_gradient_f = (Z.*d0.*dexp.*mf.^(dexp - 1))./(Z + 3) - d0.*dexp.*mf.^(dexp - 1) - (W.*a3FM.*a0W.*mf.^(a3FM - 1).*mw.^a3WM)./(1000000000.*((F.*a0W.*h0W.*mf.^a3FM.*mf.^h3FM.*mw.^a3WM.*mw.^h3WM)./1000000 + 1)) + (W.*a0W.*mf.^a3FM.*mw.^a3WM.*((F.*a3FM.*a0W.*h0W.*mf.^h3FM.*mf.^(a3FM - 1).*mw.^a3WM.*mw.^h3WM)./1000000 + (F.*a0W.*h3FM.*h0W.*mf.^a3FM.*mf.^(h3FM - 1).*mw.^a3WM.*mw.^h3WM)./1000000))./(1000000000.*((F.*a0W.*h0W.*mf.^a3FM.*mf.^h3FM.*mw.^a3WM.*mw.^h3WM)./1000000 + 1).^2) - (GGE.*Z.*a0F.*mf.^a2FM.*mz.*mz.^a2ZM)./(1000.*mf.^2.*((Z.*a0F.*h0F.*mf.^a2FM.*mf.^h2FM.*mz.^a2ZM.*mz.^h2ZM)./1000 + 1)) - (GGE.*Z.*a0F.*mf.^a2FM.*mz.*mz.^a2ZM.*((Z.*a0F.*a2FM.*h0F.*mf.^h2FM.*mf.^(a2FM - 1).*mz.^a2ZM.*mz.^h2ZM)./1000 + (Z.*a0F.*h0F.*h2FM.*mf.^a2FM.*mf.^(h2FM - 1).*mz.^a2ZM.*mz.^h2ZM)./1000))./(1000.*mf.*((Z.*a0F.*h0F.*mf.^a2FM.*mf.^h2FM.*mz.^a2ZM.*mz.^h2ZM)./1000 + 1).^2) + (GGE.*Z.*a0F.*a2FM.*mf.^(a2FM - 1).*mz.*mz.^a2ZM)./(1000.*mf.*((Z.*a0F.*h0F.*mf.^a2FM.*mf.^h2FM.*mz.^a2ZM.*mz.^h2ZM)./1000 + 1));